本發(fā)明涉及過程工程技術(shù)領(lǐng)域，尤其涉及一種復(fù)雜冶金過程模擬計算方法及系統(tǒng)。

背景技術(shù)：

目前復(fù)雜體系多相平衡計算在有機試劑分離過程設(shè)計中研究的比較多，如萃取、蒸餾等等方面。而隨著計算機技術(shù)在硬件和軟件方面的進(jìn)步發(fā)展，涉及復(fù)雜化學(xué)反應(yīng)的多相平衡計算成為可能。深入研究基于熱力學(xué)原理的多相多組分反應(yīng)體系數(shù)學(xué)模型，設(shè)計高效、魯棒性好的多相平衡計算方法已成為冶金化工領(lǐng)域的一個重要課題且已經(jīng)具有廣泛的應(yīng)用基礎(chǔ)。相關(guān)多相平衡計算的理論方法對了解工藝機理十分重要，且其計算結(jié)果是相應(yīng)過程操作、優(yōu)化、設(shè)備設(shè)計的基礎(chǔ)。

工業(yè)化生產(chǎn)過程向著連續(xù)化、大型化、復(fù)雜化的方向發(fā)展，各類復(fù)雜工程問題的優(yōu)化計算方法成為人們亟需解決的問題之一?；谧钚〖妓棺杂赡芊椒ǖ亩嘞嗥胶饽Ｐ偷那蠼饧礆w屬于工程類的優(yōu)化問題。在求解實際優(yōu)化問題過程當(dāng)中，所建立的數(shù)學(xué)模型越來越復(fù)雜，且維數(shù)越來越高，很多時候其解析性質(zhì)難以獲知，而此時傳統(tǒng)的確定性優(yōu)化算法對初值要求較高，很難實現(xiàn)對這些問題的求解。此外，局部最優(yōu)解的個數(shù)往往在問題規(guī)模增大時會迅速增加，而從大量的局部最優(yōu)解中獲得其中的全局最優(yōu)解是非常具有挑戰(zhàn)性。因而需要尋求更加高效的方法和機制實現(xiàn)優(yōu)化問題中全局最優(yōu)解的求解。

技術(shù)實現(xiàn)要素：

本發(fā)明目的在于公開一種復(fù)雜冶金過程模擬計算方法及系統(tǒng)，以快速獲得入爐元素在冶煉過程達(dá)平衡時各相中的分配。

為實現(xiàn)上述目的，本發(fā)明公開了一種復(fù)雜冶金過程模擬計算方法，包括：

以反應(yīng)體系總的吉布斯自由能函數(shù)為該數(shù)學(xué)模型的目標(biāo)函數(shù)，以輸入和輸出冶煉過程體系中各種元素的質(zhì)量相等為約束條件，建立多相平衡數(shù)學(xué)模型：

min f(x)

st.A·x＝b

x＞0

矩陣A是由各相中各組分系數(shù)組成的原子系數(shù)矩陣，對應(yīng)的b是由進(jìn)入反應(yīng)體系中各元素的總摩爾量組成的列向量，x為各相中各組分的摩爾數(shù)，f(x)為反應(yīng)體系總吉布斯自由能；

采用粒子群算法求解多相平衡下各相中組分的摩爾數(shù)；所述粒子群算法包括：

第一步：設(shè)置粒子群優(yōu)化算法的參數(shù)和要求的精度，前述粒子群優(yōu)化算法的參數(shù)包括粒子數(shù)，迭代次數(shù)上限，以及速度由以鄰域最優(yōu)值過渡到全局最優(yōu)值的迭代分界點；同時獲取約束矩陣A以及b；

第二步：獲取種群中粒子位置區(qū)間范圍；

第三步：在粒子位置區(qū)間范圍內(nèi)隨機初始化粒子群的位置信息，使得粒子的位置位于Ax＝b的超平面內(nèi)，同時把兩次隨機初始化的位置信息相減賦給對應(yīng)種群粒子的速度；

第四步：根據(jù)迭代的更新機制，在所述迭代分界點之前利用鄰域最優(yōu)值對速度進(jìn)行更新，在所述迭代分界點及之后利用全局最優(yōu)值對速度進(jìn)行更行；以及依據(jù)當(dāng)前的速度獲取更新步長用以更新粒子的位置信息，并更新種群及粒子的歷史最優(yōu)位置；

第五步：檢驗迭代次數(shù)是否達(dá)到最初給定值，如達(dá)到，種群更新結(jié)束，否則返回第四步繼續(xù)運行，得到平衡時各相中各組分的摩爾數(shù)。

本發(fā)明還公開一種執(zhí)行上述方法的相應(yīng)系統(tǒng)，至少包括：

多相平衡建模模塊，用于以反應(yīng)體系總的吉布斯自由能函數(shù)為該數(shù)學(xué)模型的目標(biāo)函數(shù)，以輸入和輸出冶煉過程體系中各種元素的質(zhì)量相等為約束條件，建立多相平衡數(shù)學(xué)模型：

min f(x)

st.A·x＝b

x＞0

矩陣A由各相中各組分系數(shù)組成的原子系數(shù)矩陣，對應(yīng)的b是由進(jìn)入反應(yīng)體系中各元素的總摩爾量組成的列向量，x為各相中各組分的摩爾數(shù)，f(x)為反應(yīng)體系總吉布斯自由能；

算法處理模塊，用于采用粒子群算法求解多相平衡下各相中組分的摩爾數(shù)；所述粒子群算法包括：

第一步：設(shè)置粒子群優(yōu)化算法的參數(shù)和要求的精度，前述粒子群優(yōu)化算法的參數(shù)包括粒子數(shù)，迭代次數(shù)上限，以及速度由以鄰域最優(yōu)值過渡到全局最優(yōu)值的迭代分界點；同時獲取約束矩陣A以及b；

第二步：獲取種群中粒子位置區(qū)間范圍；

第三步：在粒子位置區(qū)間范圍內(nèi)隨機初始化粒子群的位置信息，使得粒子的位置位于Ax＝b的超平面內(nèi)，同時把兩次隨機初始化的位置信息相減賦給對應(yīng)種群粒子的速度；

第四步：根據(jù)迭代的更新機制，在所述迭代分界點之前利用鄰域最優(yōu)值對速度進(jìn)行更新，在所述迭代分界點及之后利用全局最優(yōu)值對速度進(jìn)行更行；以及依據(jù)當(dāng)前的速度獲取更新步長用以更新粒子的位置信息，并更新種群及粒子的歷史最優(yōu)位置；

第五步：檢驗迭代次數(shù)是否達(dá)到最初給定值，如達(dá)到，種群更新結(jié)束，否則返回第四步繼續(xù)運行，得到平衡時各相中各組分的摩爾數(shù)。

本發(fā)明具有以下有益效果：

將粒子群算法用于冶金過程多相平衡組分預(yù)測，種群中粒子速度的每一維分量是同步進(jìn)行更新的，這種機制迭代的粒子能夠一直保持在可行域中進(jìn)行解的探索，可以高效實現(xiàn)優(yōu)化問題中全局最優(yōu)解的求解。而且，粒子速度的更新是分段進(jìn)行的；在迭代分界點之前，例如，前90％迭代周期內(nèi)，利用鄰域最優(yōu)值對速度進(jìn)行更新，增大了種群中粒子的多樣性；而在迭代分界點及之后，例如，后10％迭代周期內(nèi)，利用全局最優(yōu)值對速度進(jìn)行更行，可以使得種群中的粒子盡快的收斂到全局最優(yōu)解附近。藉此，可使得種群在迭代初期不易早熟，而在迭代末期能夠快速收斂至全局最優(yōu)解。

下面將參照附圖，對本發(fā)明作進(jìn)一步詳細(xì)的說明。

附圖說明

構(gòu)成本申請的一部分的附圖用來提供對本發(fā)明的進(jìn)一步理解，本發(fā)明的示意性實施例及其說明用于解釋本發(fā)明，并不構(gòu)成對本發(fā)明的不當(dāng)限定。在附圖中：

圖1是本發(fā)明實施例1公開的復(fù)雜冶金過程模擬計算方法流程圖；

圖2是本發(fā)明實施例1公開的粒子群算法的流程圖；

圖3是本發(fā)明實施例1公開的具體運算例中矩陣A的部分示意圖；

圖4是本發(fā)明實施例1公開的具體運算例中b的部分示意圖；

圖5是本發(fā)明實施例1公開的具體運算例中的鄰域拓?fù)浣Y(jié)構(gòu)圖；其中圖5(a)種群環(huán)形拓?fù)浣Y(jié)構(gòu)示意圖，圖5(b)為鄰域大小為2種群拓?fù)浣Y(jié)構(gòu)示意圖；

圖6是本發(fā)明實施例2公開的迭代次數(shù)與適應(yīng)度值的關(guān)系圖；

圖7是本發(fā)明實施例3公開的迭代次數(shù)與適應(yīng)度值的關(guān)系圖。

具體實施方式

以下結(jié)合附圖對本發(fā)明的實施例進(jìn)行詳細(xì)說明，但是本發(fā)明可以由權(quán)利要求限定和覆蓋的多種不同方式實施。

實施例1

本實施例以銅冶煉為例，公開一種復(fù)雜冶金過程模擬計算方法，如圖1所示，包括：

步驟S1、以反應(yīng)體系總的吉布斯自由能函數(shù)為該數(shù)學(xué)模型的目標(biāo)函數(shù)，以輸入和輸出冶煉過程體系中各種元素的質(zhì)量相等為約束條件，建立多相平衡數(shù)學(xué)模型。

在銅冶煉工藝過程中，其反應(yīng)體系的總吉布斯自由能可由下式(1)和(2)計算：

G(n,T,P)＝G₁(n,T,P)_Matte+G₂(n,T,P)_Slag+G₃(n,T,P)_Gas (1)

上式(1)中，G(n,T,P)，G₁(n,T,P)_Matte，G₂(n,T,P)_Slag，G₃(n,T,P)_Gas分別為反應(yīng)體系中總吉布斯自由能，以及銅锍相，爐渣相和氣相吉布斯自由能。而式(2)中，N_p和N_c分別表示體系總相數(shù)和各相中組分?jǐn)?shù)，n_ij為j相中組分i的摩爾量，μ_ij表示在反應(yīng)體系條件下j相中組分i的化學(xué)勢。ΔG_ij^o表示在系統(tǒng)溫度、標(biāo)準(zhǔn)壓強下j相中組分i的吉布斯生成自由能，R是理想氣體常數(shù)，T表示開爾文溫度，f_ij是在j相中組分i的偏逸度，是j相中組分i在參考狀態(tài)下的逸度。

受質(zhì)量守恒和組分摩爾量為非負(fù)的約束，這種約束條件在冶金化工領(lǐng)域常見。在銅火法冶煉工藝過程中，由于CaO、MgO、Al₂O₃、SiO₂這類惰性組分直接全部進(jìn)入渣相，而不參與反應(yīng)，因而爐渣中這部分組分的摩爾量可以由其進(jìn)入量而直接得到。而混合礦中所帶入的水分，全部轉(zhuǎn)化為水蒸氣而進(jìn)入氣相中，因而氣相中的水蒸氣摩爾量也可以直接得到。同時氣相中的氮氣呈惰性物質(zhì)，可以從進(jìn)入反應(yīng)體系中的富氧計算得到。其他的元素如Cu、Fe、S、O、As、Sb、Bi、Pb、Zn等由于平衡時候在三相中均有對應(yīng)的組分存在，因而其質(zhì)量守恒條件滿足，進(jìn)入體系的元素總量等于其在三相中分配元素的總量。其質(zhì)量守恒表達(dá)式可以用下式(3)表示：

Ax＝b (3)

其中A為由對應(yīng)組分化學(xué)式組成的原子矩陣；x為達(dá)到平衡時，由三相中各組分的摩爾量而組成的列向量，即為最終需要求解的變量；而b為由進(jìn)入體系中元素的總量而組成的列向量。同時由于x中各組分的摩爾量不允許出現(xiàn)負(fù)數(shù)，所以必須同時保證x>0。

藉此，本實施例所建立多相平衡數(shù)學(xué)模型為：

min f(x)

st.A·x＝b

x＞0。

步驟S2、采用粒子群算法求解多相平衡下各相中組分的摩爾數(shù)。

如圖2所示，該步驟可細(xì)分為：

步驟S21、設(shè)置粒子群優(yōu)化算法的參數(shù)和要求的精度，前述粒子群優(yōu)化算法的參數(shù)包括粒子數(shù)，迭代次數(shù)上限，以及速度由以鄰域最優(yōu)值過渡到全局最優(yōu)值的迭代分界點；同時獲取約束矩陣A以及b。

步驟S22、獲取種群中粒子位置區(qū)間范圍。換言之，即確定變量x中第j維分量的最大上限值，可以通過下式確定：

x_jmax＝min(b_i/A_ij|A_ij≠0) (4)

其中，b_i為體系中相應(yīng)元素的含量，A_ij為矩陣A中相應(yīng)的原子系數(shù)。

步驟S23、在粒子位置區(qū)間范圍內(nèi)隨機初始化粒子群的位置信息，使得粒子的位置位于Ax＝b的超平面內(nèi)，同時把兩次隨機初始化的位置信息相減賦給對應(yīng)種群粒子的速度。

可選的，初始化位置和速度信息包括：

從矩陣A的n列中選出線性無關(guān)的m列，用B表示這m階方陣，用C表示A中剩下的(n-m)列而組成的m×(n-m)子矩陣，對應(yīng)的變量x可以相應(yīng)的分解為x＝[x_B；x_C]，則可以轉(zhuǎn)化為：

A·x＝[B,C]·[x_B；x_C]＝B·x_B+C·x_C＝b (5)

則x_B即可通過下式得到：

x_B＝B^-1·b-B^-1·C·x_C (6)

通過上式中隨機賦值非基變量x_C于區(qū)間[0,x_Cmax]中，對于求解出的基變量x_B，檢驗其中x_B的每一維分量是否均大于0，如果是，則x＝[x_B；x_C]即成功初始化在約束超平面Ax＝b內(nèi)；如果x_B中每一維分量并不都大于0，則需要重新隨機賦值x_C于區(qū)間[0,x_Cmax]中直至求解出的x_B中每一維分量都大于0；其中，x_Cmax為粒子變量各相應(yīng)維度的最大上限值。

步驟S24、根據(jù)迭代的更新機制，在迭代分界點之前利用鄰域最優(yōu)值對速度進(jìn)行更新，在迭代分界點及之后利用全局最優(yōu)值對速度進(jìn)行更行；以及依據(jù)當(dāng)前的速度獲取更新步長用以更新粒子的位置信息，并更新種群及粒子的歷史最優(yōu)位置。

優(yōu)選的，在迭代更新時，還可以對種群中的粒子增加位置和速度擾動。例如：

位置擾動是通過三個粒子位置向量線性組合生成一個臨時的粒子位置pbest_temp實現(xiàn)的，計算公式如下：

pbest_temp＝pbest_i+r·(pbset_rand1-pbset_rand2) (7)

其中兩個粒子位置向量pbest_rand1和pbest_rand2是從當(dāng)前粒子最優(yōu)位置池中隨機選擇，r是區(qū)間(0,1)之間的隨機數(shù)，擾動位置pbest_temp會同相應(yīng)的當(dāng)前位置pbest_i進(jìn)行比較，如果擾動位置pbest_temp相比當(dāng)前位置pbest_i在種群中有更好的適應(yīng)度值，則粒子當(dāng)前的最優(yōu)位置pbest_i就會被pbest_temp替換。

速度擾動是當(dāng)前速度通過一個稱為速度狀態(tài)矩陣進(jìn)行線性變化生成一個臨時的粒子速度v_temp實現(xiàn)的，如下式所示：

v_temp＝v_i·v/||v|| (8)

速度狀態(tài)矩陣v中的行是由從當(dāng)前粒子速度池中隨機選取出的m個粒子速度向量組成；在臨時速度基礎(chǔ)上，更新當(dāng)前粒子的位置產(chǎn)生一個新的pbest_temp如下式：

pbest_temp＝pbest_i+λ·v_temp (9)

最后，比較當(dāng)前粒子的最優(yōu)位置與新產(chǎn)生的pbest_temp在種群中的適應(yīng)度值，如果pbest_temp較優(yōu)，則用pbest_temp替代原先的最優(yōu)位置pbest_i；其中λ為步長。

其中，在該步驟S24所進(jìn)行種群中的粒子i更新位置的最大步長λ_{i max}，通過下式獲取：

λ_imax＝min(x_jmax-x_ij/v_ij|v_ij＞ε,-x_ij/v_ij|v_ij＜-ε,0|-ε＜v_ij＜ε) (10)

其中ε為計算精度，值為用來避免0作為除數(shù)但約等于0的正數(shù)，x_jmax為變量x中第j維分量的最大上限值，x_ij為粒子i的j維分量的位置，v_ij為粒子i的j維分量的速度。

換言之：上式(9)中，i表示第i個粒子，j表示該粒子的第j維。ε為計算精度用來避免0作為除數(shù)，其值約等于0的正數(shù)，其中：

當(dāng)v>ε(其含義可視為v>0)時，視為粒子向正方向運動，其邊界為x_i(t+1)＝x_max；

當(dāng)v<-ε(其含義可視為v<0)時，視為粒子向負(fù)方向運動，其邊界為x_i(t+1)＝0；

當(dāng)ε<v<-ε(其含義可視為v＝0)時，視為粒子不運動，此時不需對粒子步長進(jìn)行約束。

優(yōu)選的，本實施例中，種群中粒子i的步長系數(shù)λ_i按如下公式進(jìn)行更新：

在上述迭代更新步驟中，所謂“適應(yīng)度”是將求解出的粒子代入反應(yīng)體系總吉布斯自由f(x)中，評估標(biāo)準(zhǔn)為：f(x)越小，其適應(yīng)度越好。

步驟S25、檢驗迭代次數(shù)是否達(dá)到最初給定值，如達(dá)到，種群更新結(jié)束，否則返回步驟S24繼續(xù)運行，得到平衡時各相中各組分的摩爾數(shù)。

為便于本領(lǐng)域技術(shù)人員對上述過程做進(jìn)一步理解，本發(fā)明以具體的運算例對上述過程進(jìn)行復(fù)述：

在銅冶煉工程中，x為Cu₂S、Cu、FeS、FeO、Fe₃O₄、Pb、PbS、ZnS、As、Sb、Bi，F(xiàn)eO、Cu₂O、Fe₃O₄、FeS、PbO、ZnO、As₂O₃、Sb₂O₃、Bi₂O₃、SiO₂、CaO、MgO、Al₂O₃，SO₂、S₂、O₂、N₂、H₂O、PbO、PbS、Zn、ZnS、As₂、AsO、AsS、SbO、SbS、BiO等這些物質(zhì)在熔煉過程中的含量。其中，銅锍相中主要包括Cu₂S、Cu、FeS、FeO、Fe₃O₄、Pb、PbS、ZnS、As、Sb、Bi等組分；渣相中主要包括FeO、Cu₂S、Cu₂O、Fe₃O₄、FeS、PbO、ZnO、As₂O₃、Sb₂O₃、Bi₂O₃、SiO₂、CaO、MgO、Al₂O₃等組分；氣相中主要包括SO₂、S₂、O₂、N₂、H₂O、PbO、PbS、Zn、ZnS、As₂、AsO、AsS、SbO、SbS、BiO等組分。

系數(shù)矩陣A為上述物質(zhì)組成的系數(shù)矩陣，截圖的部分矩陣A如圖3所示，其中，矩陣的行是某一元素在上述物質(zhì)中的原子系數(shù)(即腳標(biāo))，以金屬Cu為例，它在A中表示為：[2 1 2]。

如圖4所示，b為體系中某元素的含量(物質(zhì)的量)，可根據(jù)投入物料計算得到；以金屬Cu為例，在加料量66噸，銅含量24.352％時，體系中Cu：66*1000*1000/64＝253107.40157480314mol。

由于f(x)函數(shù)對應(yīng)的等式方程為公式(2)，而公式(2)及公式(3)組成的多相平衡數(shù)學(xué)模型是欠定方程組，理論上有很多解，這里就是利用粒子群算法求解該方程，使其解滿足方程(3)及x>0的同時使公式(2)達(dá)到最小值。粒子的位置即為方程的解，也即為體系中物質(zhì)的含量。

由于所有元素組成的Ax＝b太復(fù)雜，這里選取元素銅、鐵為例說明：

體系中含銅的物質(zhì)有三種：Cu₂S、Cu、Cu₂O，設(shè)其含量為x₁、x₂、x₃，(其他含銅化合物含量少，忽略不計)

體系中含鐵的物質(zhì)有三種：FeS、FeO、Fe₃O₄,設(shè)其含量為x₄、x₅、x₆,(其他含鐵化合物含量少，忽略不計)

其Ax＝b可以表示為：

將上述方程(12)展開，得：

對于公式(a)，要想使所有的x均大于0，則每個x必有一個最大值，該最大值必須小于常數(shù)項與系數(shù)的比值，以x₁為例，其最大值為253107.40/2，因為如果x₁大于該值，其他x必須取負(fù)值才能滿足上述等式方程。

上述是比較特殊的情況，矩陣A中每一列只有一個值，若矩陣A中每一列有兩個或更多的值，如下：

將上述方程展開，得

在該方程組內(nèi)，如按上述規(guī)則取x₁的最大值，會對應(yīng)兩個值(253107.40/2,315385.71/5)，則x₁取最小的一個(315385.71/5)作為最大值，才能使所有的x取值非負(fù)。

x_1max＝min{253107.40/2,315385.71/5}＝315385.71/5 (16)

上述式(15)及式(16)可用于解釋上述步驟S22中的獲取種群中粒子位置區(qū)間范圍。下面對粒子群的初始化做進(jìn)一步描述：

為保證計算的速度和精度，在計算過程中會使用若干粒子(數(shù)量與實際問題的復(fù)雜程度有關(guān))組成的粒子群，由于每個粒子的初始化和迭代過程類似，這里以一個粒子(假設(shè)問題只考慮Fe、Cu)為例，該粒子共有6維，每一維即為一個物種(x_i)。

從上式(13)中的矩陣A中取出m(m＝2)線性無關(guān)列，構(gòu)成基向量剩余的構(gòu)成非基向量從新調(diào)整矩陣使其仍滿足

由線性代數(shù)知識，在上述范圍內(nèi)任意給定x_C的值，可以求得:

x_B＝B^-1·b-B^-1·C·x_C (17)

其中，x_Cmax為粒子變量各相應(yīng)維度的最大上限值。通過上式(17)中隨機賦值非基變量x_C于區(qū)間[0,x_Cmax]中，對于求解出的基變量x_B，檢驗其中x_B的每一維分量是否均大于0，如果是，則x＝[x_B；x_C]即成功初始化在約束超平面Ax＝b內(nèi)；如果x_B中每一維分量并不都大于0，則需要重新隨機賦值x_C于區(qū)間[0,x_Cmax]中直至求解出的x_B中每一維分量都大于0；即得到了問題的初始解(一個粒子的初始位置)。粒子在運動過程中，兩次位置之差即為粒子運動的速度v，因此，可以對每個粒子初始化兩次，將兩次位置之差記為粒子的初始化速度。

v＝x(t+1)-x(t) (18)

由于上述初始化滿足Ax＝b，即為方程組的解，但是并不一定能滿足目標(biāo)函數(shù)值達(dá)到最小，因此需要對粒子群的位置進(jìn)行更新，使其滿足Ax＝b，且x>0的前提下，使目標(biāo)函數(shù)值f(x)最小。為此，關(guān)于上述步驟S24中的例子速度及位置更新以及擾動等，通過下述(一)至(四)詳述如下：

(一)、種群迭代步長求解

為了避免種群粒子在下次速度更新完之后越過可行解搜索域，需要對初始速度進(jìn)行相應(yīng)的速度大小限制，通過增加步長系數(shù)λ即可滿足粒子每次速度更新后不會越過可行解搜索域。種群中粒子i的最大步長系數(shù)λ_{i max}可以通過上式(10)得到：

λ_imax＝min(x_jmax-x_ij/v_ij|v_ij＞ε,-x_ij/v_ij|v_ij＜-ε,0|-ε＜v_ij＜ε) (10)

如果所有的粒子在更新的時候都使用λ_{i max}，相應(yīng)的實驗表明種群的局部搜索性能會逐漸減小，因而在本實施例，當(dāng)粒子的λ_{i max}>1時候，這部分粒子的更新步長系數(shù)重新被置為1，而其他粒子保持不變。因而種群中粒子i的步長系數(shù)λ_i按如下公式(11)進(jìn)行更新。

(二)、種群拓?fù)浣Y(jié)構(gòu)設(shè)計

實際過程中，對于復(fù)雜難優(yōu)化的問題，要求種群有好的全局搜索能力，而當(dāng)種群定位到最優(yōu)解附近位置時，要求加強局部探索能力，即獲取精度更高的最優(yōu)解。好的種群拓?fù)浣Y(jié)構(gòu)能夠有效的保證搜索初期在全局范圍內(nèi)避免種群早熟的探尋能力，以及搜索后期在局部范圍內(nèi)更精確的探尋能力。

為此，如附圖5(a)所示，種群拓?fù)浣Y(jié)構(gòu)可采用環(huán)形結(jié)構(gòu)；本實施例設(shè)計了一種名為單鏈環(huán)形鄰域拓?fù)浣Y(jié)構(gòu)的粒子群，一種典型的鄰域大小為2(鄰域大小為所考慮的周圍粒子的數(shù)量，不代表距離)的種群拓?fù)浣Y(jié)構(gòu)如附圖5(b)所示，對于種群中粒子k，鄰域由種群中粒子編號為k+1，k-2，k+3，k-4…k+n-1和k-n的N個粒子組成；這種拓?fù)浣Y(jié)構(gòu)，能夠有效避免編號相鄰粒子間的相互吸引。通過實施上述選擇鄰域的策略，種群中粒子i的鄰域最佳粒子記為lbest_i。

(三)、種群全局搜索與局部探索能力平衡機制

從全局的視角出發(fā)平衡種群的全局搜索能力和局部探索能力，對種群的權(quán)重因子影響和速度更新策略進(jìn)行了充分的考慮。一般，迭代初期需要保持種群有較強的全局搜索能力，種群的粒子能夠盡可能分散在搜索域的各個角落；而迭代后期需要保持種群有較強的局部探索能力，種群的粒子能夠盡快收斂到全局最優(yōu)值附近，滿足探索更高精度的優(yōu)化解。

粒子群算法中，高的權(quán)重因子就意味著粒子能夠快速的在搜索空間內(nèi)移動，保證粒子能夠達(dá)到搜索空間的任意位置。因而權(quán)重因子在迭代初期應(yīng)該保持高的水平，而隨著迭代次數(shù)的增加，權(quán)重因子w逐漸降低到零。在本發(fā)明中的粒子群優(yōu)化算法中，權(quán)重因子的更新遵循下面的公式(19)：

w＝exp(-30·(k/MaxIter)¹⁰) (19)

在上述公式中，MaxIter是種群的總迭代次數(shù)，k是當(dāng)前的迭代次數(shù)。藉此，權(quán)重因子隨著迭代次數(shù)的變化呈現(xiàn)一條S曲線進(jìn)行變化。

在本實施例中，粒子速度的更新是分段進(jìn)行的。例如：在前90％迭代周期內(nèi)，種群中粒子的速度更新遵循公式(20)，為了增大種群中粒子的多樣性；而在后10％迭代周期內(nèi)，種群中粒子的速度更新遵循公式(21)，為了使得種群中的粒子盡快的收斂到最優(yōu)解附近。因此上述的策略能夠使得種群在迭代初期不易早熟，而在迭代末期能夠快速收斂至全局最優(yōu)解。種群中粒子速度的每一維分量是同步進(jìn)行更新的，這種機制迭代的粒子能夠一直保持在可行域中進(jìn)行解的探索。

v_i(t+1)＝w·v_i(t)+r₁·[pbest_i(t)-x_i(t)]+r₂·[lbest(t)-x_i(t)] (20)

v_i(t+1)＝w·v_i(t)+r₁·[pbest_i(t)-x_i(t)]+r₂·[gbest(t)-x_i(t)] (21)

其中，i表示第i個粒子，(t)表示第t次更新，(t+1)表示第t+1次更新，x表示粒子位置，v表示粒子速度。pbest_i為粒子當(dāng)前的最優(yōu)位置，lbest_i為粒子鄰域內(nèi)的最優(yōu)位置，gbest_i為粒子群內(nèi)當(dāng)前的最優(yōu)位置，r₁和r₂是區(qū)間(0,1)之間兩相互獨立的隨機數(shù)。

粒子的位置更新可以用下式(21)表示：

x_i(t+1)＝x_i(t)+λ_i·v_i(t+1) (22)

(四)、種群粒子位置和速度擾動

對于當(dāng)前粒子i最優(yōu)位置pbest_i以及速度v_i添加擾動能夠在保持種群多樣性的同時，引導(dǎo)粒子飛向更優(yōu)的位置而不破壞種群自身的組織結(jié)構(gòu)。

種群中粒子位置擾動是通過三個粒子位置向量線性組合生成一個臨時的粒子位置pbest_temp實現(xiàn)的，如式(7)所示。其中兩個粒子位置向量是從的當(dāng)前粒子最優(yōu)位置池中隨機選擇，r是區(qū)間(0,1)之間的隨機數(shù)。通過對當(dāng)前位置pbest_i添加擾動，使得產(chǎn)生的pbest_temp仍然滿足A·pbest_temp＝b。之后擾動位置pbest_temp會同相應(yīng)的當(dāng)前位置pbest_i進(jìn)行比較，如果擾動位置pbest_temp相比當(dāng)前位置pbest_i在種群中有更好的適應(yīng)度值，則粒子當(dāng)前的最優(yōu)位置pbest_i就會被pbest_temp替換。

pbest_temp＝pbest_i+r·(pbset_rand1-pbset_rand2) (7)

種群中對粒子速度的擾動是把當(dāng)前速度通過一個稱為速度狀態(tài)矩陣進(jìn)行線性變化生成一個臨時的粒子速度v_temp實現(xiàn)的，如式8)所示。

v_temp＝v_i·v/||v|| (8)

矩陣v即為速度狀態(tài)矩陣，是一個m×m的方陣，m即為粒子群的大小。速度狀態(tài)矩陣v中的行是由從當(dāng)前粒子速度池中隨機選取出的m個粒子速度向量組成。因而當(dāng)前的粒子速度v_i在速度狀態(tài)矩陣v轉(zhuǎn)化下，能夠向任何方向進(jìn)行偏移，種群的局部搜索能力能夠得到大大的提高。而經(jīng)過此擾動，生成的v_temp仍然能夠滿足A·v_temp＝0條件。在此基礎(chǔ)上，根據(jù)式(10)和(11)確定每次迭代的步長，更新當(dāng)前粒子的位置產(chǎn)生一個新的pbest_temp如式(9)所示。最后，比較當(dāng)前粒子的最優(yōu)位置與新產(chǎn)生的pbest_temp在種群中的適應(yīng)度值，如果pbest_temp較優(yōu)，則用pbest_temp替代原先的最優(yōu)位置pbest_i。

pbest_temp＝pbest_i+λ·v_temp (9)

綜上，本發(fā)明適用于求解帶超高維數(shù)線性約束，且非凸目標(biāo)函數(shù)優(yōu)化問題的改進(jìn)粒子群優(yōu)化算法(HLPSO)的主旨如下：最先需要獲取原子矩陣A和列向量b，以及算法的相關(guān)設(shè)置參數(shù)。由于約束條件變量x非負(fù)的要求，同時A和b中的成員均為非負(fù)，對于變量x來說自然存在上界。保持變量x不超出上界保證了后續(xù)的迭代操作均約束在Ax＝b超平面內(nèi)。之后即為最重要的兩個步驟：種群初始化和種群過程迭代。

HLPSO算法中種群的初始化包括粒子位置初始化和速度初始化。由于粒子初始位置是隨機生成的，因而有一定的概率初始化失敗，直至初始化所有粒子的位置x在0與上界之間的區(qū)間。而粒子的速度初始化即為兩次粒子隨機初始化的位置之差。在所有粒子均初始化成功后，選擇種群中所有粒子位置中的最優(yōu)位置為全局最優(yōu)位置。上述種群初始化成功后，即進(jìn)入種群過程迭代。

種群迭代過程即粒子學(xué)習(xí)過程，逐步靠近全局最優(yōu)值的過程。每次迭代過程中，獲取種群全局及粒子鄰域最優(yōu)位置用于粒子速度更新，而速度更新過程中需要對速度做相應(yīng)的限制，以免在位置更新時候躍出可行搜索區(qū)域。之后更新種群及粒子的歷史最優(yōu)位置，在此基礎(chǔ)上從粒子池中隨機選擇粒子數(shù)據(jù)添加對當(dāng)前粒子的位置和速度的擾動，以避免迭代前期陷入局部最優(yōu)值，同時在迭代后期能夠加強粒子局部探索能力。

迭代次數(shù)達(dá)到后，種群更新結(jié)束，而最后種群中的歷史全局最優(yōu)位置極可能為所求解問題的全局最優(yōu)值。

實施例2：

為了測試算法的性能，從文獻(xiàn)中選取了如下式(23)的測試函數(shù)，其函數(shù)形式同本設(shè)計銅冶煉過程多相平衡數(shù)學(xué)模型類型相類似，具有一定的代表性。

其中c₁＝-6.089,c₂＝-17.164,c₃＝-34.054,c₄＝-5.914,c₅＝-24.721,c₆＝-14.896,c₇＝-24.100,c₈＝-10.708,c₉＝-26.662,c₁₀＝-22.179。

COPSO算法是由Aguirre A H，Zavala A M et al.設(shè)計提出的，ISRES算法是由Runarsson和Yao提出的。其與HLPSO針對上述測試函數(shù)的測試結(jié)果比較如下表1所示。

表1：

文獻(xiàn)中給出的精確解為：

x*＝(0.0407，0.1477，0.7832，0.0013，0.4853，0.0007，0.0274，0.018，0.0373，0.0968)

而使用HLPSO每次都能收斂到上述最優(yōu)解附近，十次測試的平均最優(yōu)解為：x＝(0.0407，0.1477，0.7831，0.0014，0.4853，6.959e-04，0.0274，0.018，0.0373，0.0968)，由此可見與文獻(xiàn)中給出的精確解幾乎一樣，說明HLPSO算法的是可行有效的。

由上述比較可知，HLPSO算法在處理約束方面的內(nèi)在特性，種群能夠始終保持在可行域內(nèi)，因而搜索效率相比于COPSO和ISRES性能大大提高，其穩(wěn)定性也好于其他兩種算法，由圖6可以知道，HLPSO算法很快的就能夠搜索到全局最優(yōu)值附近，在總計100次迭代內(nèi)，大約40次迭代后就收斂到了全局最優(yōu)值，收斂速度塊，具有較高的搜索效率。因而HLPSO算法是適合針對本發(fā)明的所提出的銅冶煉過程多相平衡數(shù)學(xué)模型的求解。

實施例3：

在國內(nèi)某銅冶煉廠，采集了相關(guān)生產(chǎn)實踐過程中的基礎(chǔ)參數(shù)，建立多相平衡數(shù)學(xué)模型。

通過生產(chǎn)調(diào)研，冰銅和爐渣相的溫度約1200℃，而氣相溫度為1300℃；銅锍相中主要包括Cu₂S、Cu、FeS、FeO、Fe₃O₄、Pb、PbS、ZnS、As、Sb、Bi等組分；渣相中主要包括FeO、Cu₂S、Cu₂O、Fe₃O₄、FeS、PbO、ZnO、As₂O₃、Sb₂O₃、Bi₂O₃、SiO₂、CaO、MgO、Al₂O₃等組分；氣相中主要包括SO₂、S₂、O₂、N₂、H₂O、PbO、PbS、Zn、ZnS、As2、AsO、AsS、SbO、SbS、BiO等組分。

在上述基礎(chǔ)上通過熱力學(xué)多相平衡理論建立數(shù)學(xué)模型，確定了優(yōu)化目標(biāo)函數(shù)以及約束條件，并對其進(jìn)行計算求取最終平衡時候各相中的組分含量。

從國內(nèi)某以銅冶煉廠獲取了其在2014年穩(wěn)定工況下某段時間內(nèi)的入爐混合銅精礦成分表及其對應(yīng)的操作工藝條件。其中表2為入爐銅精礦的成分表，而工藝操作過程參數(shù)如表3所示。

表2：

續(xù)

表3：

在上述生產(chǎn)工況下，選取處理得到的工業(yè)生產(chǎn)數(shù)據(jù)同模擬計算的結(jié)果進(jìn)行了對比，如下表4和表5所示。

表4：

表5：

通過以上模擬計算數(shù)據(jù)和工業(yè)生產(chǎn)數(shù)據(jù)對比分析結(jié)果可知，所建立的銅冶煉過程多相平衡模型是可靠的，可用于仿真預(yù)測銅冶煉工藝過程中多元素的分配行為以及各組分之間的相互關(guān)系。

設(shè)置HLPSO算法中種群大小為200，而迭代次數(shù)為1000，重復(fù)進(jìn)行相同的測試10次。10次重復(fù)測試實驗的結(jié)果如下表6所示。

表6：

由上述測試可知，針對銅冶煉過程多相平衡數(shù)學(xué)模型，HLPSO算法能夠勝任上述模型的計算，計算成功率為90％，幾乎每次都能收斂至全局最優(yōu)值附近。圖7說明了HLPSO算法很快就能夠搜索到全局最優(yōu)值附近，在總計1000次迭代內(nèi)，大約200次迭代后就收斂到了全局最優(yōu)值，收斂速度塊，具有較高的搜索效率。因而HLPSO算法能夠適用于本發(fā)明所提出的銅冶煉過程多相平衡數(shù)學(xué)模型的求解，且HLPSO算法是可靠的。

實施例4

與上述實施例相對應(yīng)的，本發(fā)明實施例還公開一種復(fù)雜冶金過程模擬計算系統(tǒng)，至少包括：

多相平衡建模模塊，用于以反應(yīng)體系總的吉布斯自由能函數(shù)為該數(shù)學(xué)模型的目標(biāo)函數(shù)，以輸入和輸出冶煉過程體系中各種元素的質(zhì)量相等為約束條件，建立多相平衡數(shù)學(xué)模型。

如上述實施例1，該多相平衡建模模塊所建立的多相平衡數(shù)學(xué)模型具體為：

min f(x)

st.A·x＝b

x＞0

上述模型中，矩陣A是由各相中各組分系數(shù)組成的原子系數(shù)矩陣，對應(yīng)的b是由進(jìn)入反應(yīng)體系中各元素的總摩爾量組成的列向量，x為各相中各組分的摩爾數(shù)，f(x)為反應(yīng)體系總吉布斯自由能；

算法處理模塊，用于采用粒子群算法求解多相平衡下各相中組分的摩爾數(shù)。所述粒子群算法包括：

第一步：設(shè)置粒子群優(yōu)化算法的參數(shù)和要求的精度，前述粒子群優(yōu)化算法的參數(shù)包括粒子數(shù)，迭代次數(shù)上限，以及速度由以鄰域最優(yōu)值過渡到全局最優(yōu)值的迭代分界點；同時獲取約束矩陣A以及b；

第二步：獲取種群中粒子位置區(qū)間范圍；

第三步：在粒子位置區(qū)間范圍內(nèi)隨機初始化粒子群的位置信息，使得粒子的位置位于Ax＝b的超平面內(nèi)，同時把兩次隨機初始化的位置信息相減賦給對應(yīng)種群粒子的速度；

第四步：根據(jù)迭代的更新機制，在所述迭代分界點之前利用鄰域最優(yōu)值對速度進(jìn)行更新，在所述迭代分界點及之后利用全局最優(yōu)值對速度進(jìn)行更行；以及依據(jù)當(dāng)前的速度獲取更新步長用以更新粒子的位置信息，并更新種群及粒子的歷史最優(yōu)位置；

第五步：檢驗迭代次數(shù)是否達(dá)到最初給定值，如達(dá)到，種群更新結(jié)束，否則返回第四步繼續(xù)運行，得到平衡時各相中各組分的摩爾數(shù)。

其中，本實施例中算法處理模塊的具體計算請參照上述實施例1，在此不做贅述。綜上，本實施例所公開的復(fù)雜冶金過程模擬計算方法及系統(tǒng)，可以高效實現(xiàn)冶金過程多相平衡組分預(yù)測等優(yōu)化問題中全局最優(yōu)解的求解。具體的：將粒子群算法用于冶金過程多相平衡組分預(yù)測，種群中粒子速度的每一維分量是同步進(jìn)行更新的，這種機制迭代的粒子能夠一直保持在可行域中進(jìn)行解的探索，可以高效實現(xiàn)優(yōu)化問題中全局最優(yōu)解的求解。而且，粒子速度的更新是分段進(jìn)行的；在迭代分界點之前，例如，前90％迭代周期內(nèi)，利用鄰域最優(yōu)值對速度進(jìn)行更新，增大了種群中粒子的多樣性；而在迭代分界點及之后，例如，后10％迭代周期內(nèi)，利用全局最優(yōu)值對速度進(jìn)行更行，可以使得種群中的粒子盡快的收斂到全局最優(yōu)解附近。藉此，可使得種群在迭代初期不易早熟，而在迭代末期能夠快速收斂至全局最優(yōu)解。

此外，優(yōu)選的，當(dāng)本實施例的上述方法及系統(tǒng)應(yīng)用于冶煉銅時，還可以進(jìn)一步將多相平衡數(shù)學(xué)模型與機械夾雜方程結(jié)合起來，對主要成分進(jìn)行機械夾雜的修正；其中：

銅锍在爐渣中的夾雜率計算公式為：

爐渣在銅锍相中的夾雜率計算公式為：

銅锍相中夾雜爐渣量的計算公式為：

爐渣相中夾雜銅锍量的計算公式為：

其中，M_slag和M_matte分別表示理論計算平衡時爐渣相和銅锍相的總質(zhì)量，M_slag和M_matte分別是進(jìn)入銅锍相中的爐渣質(zhì)量和進(jìn)入爐渣相中的銅锍質(zhì)量，和分別是爐渣在銅锍相中的夾雜系數(shù)以及銅锍在爐渣相中的夾雜系數(shù)。

以上所述僅為本發(fā)明的優(yōu)選實施例而已，并不用于限制本發(fā)明，對于本領(lǐng)域的技術(shù)人員來說，本發(fā)明可以有各種更改和變化。凡在本發(fā)明的精神和原則之內(nèi)，所作的任何修改、等同替換、改進(jìn)等，均應(yīng)包含在本發(fā)明的保護(hù)范圍之內(nèi)。