CN106611104B - 复杂冶金过程模拟计算方法及系统 - Google Patents

Abstract

本发明涉及过程工程技术领域，公开一种复杂冶金过程模拟计算方法及系统，以快速获得入炉元素在冶炼过程达平衡时各相中的分配。本发明公开的方法包括：以反应体系总的吉布斯自由能函数为该数学模型的目标函数，以输入和输出冶炼过程体系中各种元素的质量相等为约束条件，建立多相平衡数学模型，并结合机械夹杂方程对该数学模型进行修正，然后采用粒子群算法求解多相平衡下各相中组分的摩尔数；所述粒子群算法根据迭代的更新机制，在迭代分界点之前利用邻域最优值对速度进行更新，在迭代分界点及之后利用全局最优值对速度进行更行；以及依据当前的速度获取更新步长用以更新粒子的位置信息，并更新种群及粒子的历史最优位置。

Description

复杂冶金过程模拟计算方法及系统

技术领域

本发明涉及过程工程技术领域，尤其涉及一种复杂冶金过程模拟计算方法及系统。

背景技术

目前复杂体系多相平衡计算在有机试剂分离过程设计中研究的比较多，如萃取、蒸馏等等方面。而随着计算机技术在硬件和软件方面的进步发展，涉及复杂化学反应的多相平衡计算成为可能。深入研究基于热力学原理的多相多组分反应体系数学模型，设计高效、鲁棒性好的多相平衡计算方法已成为冶金化工领域的一个重要课题且已经具有广泛的应用基础。相关多相平衡计算的理论方法对了解工艺机理十分重要，且其计算结果是相应过程操作、优化、设备设计的基础。

工业化生产过程向着连续化、大型化、复杂化的方向发展，各类复杂工程问题的优化计算方法成为人们亟需解决的问题之一。基于最小吉布斯自由能方法的多相平衡模型的求解即归属于工程类的优化问题。在求解实际优化问题过程当中，所建立的数学模型越来越复杂，且维数越来越高，很多时候其解析性质难以获知，而此时传统的确定性优化算法对初值要求较高，很难实现对这些问题的求解。此外，局部最优解的个数往往在问题规模增大时会迅速增加，而从大量的局部最优解中获得其中的全局最优解是非常具有挑战性。因而需要寻求更加高效的方法和机制实现优化问题中全局最优解的求解。

发明内容

本发明目的在于公开一种复杂冶金过程模拟计算方法及系统，以快速获得入炉元素在冶炼过程达平衡时各相中的分配。

为实现上述目的，本发明公开了一种复杂冶金过程模拟计算方法，包括：

以反应体系总的吉布斯自由能函数为该数学模型的目标函数，以输入和输出冶炼过程体系中各种元素的质量相等为约束条件，建立多相平衡数学模型：

min f(x)

st.A·x＝b

x＞0

矩阵A是由各相中各组分系数组成的原子系数矩阵，对应的b是由进入反应体系中各元素的总摩尔量组成的列向量，x为各相中各组分的摩尔数，f(x)为反应体系总吉布斯自由能；

采用粒子群算法求解多相平衡下各相中组分的摩尔数；所述粒子群算法包括：

第一步：设置粒子群优化算法的参数和要求的精度，前述粒子群优化算法的参数包括粒子数，迭代次数上限，以及速度由以邻域最优值过渡到全局最优值的迭代分界点；同时获取约束矩阵A以及b；

第二步：获取种群中粒子位置区间范围；

第三步：在粒子位置区间范围内随机初始化粒子群的位置信息，使得粒子的位置位于Ax＝b的超平面内，同时把两次随机初始化的位置信息相减赋给对应种群粒子的速度；

第四步：根据迭代的更新机制，在所述迭代分界点之前利用邻域最优值对速度进行更新，在所述迭代分界点及之后利用全局最优值对速度进行更行；以及依据当前的速度获取更新步长用以更新粒子的位置信息，并更新种群及粒子的历史最优位置；

第五步：检验迭代次数是否达到最初给定值，如达到，种群更新结束，否则返回第四步继续运行，得到平衡时各相中各组分的摩尔数。

本发明还公开一种执行上述方法的相应系统，至少包括：

多相平衡建模模块，用于以反应体系总的吉布斯自由能函数为该数学模型的目标函数，以输入和输出冶炼过程体系中各种元素的质量相等为约束条件，建立多相平衡数学模型：

min f(x)

st.A·x＝b

x＞0

矩阵A由各相中各组分系数组成的原子系数矩阵，对应的b是由进入反应体系中各元素的总摩尔量组成的列向量，x为各相中各组分的摩尔数，f(x)为反应体系总吉布斯自由能；

算法处理模块，用于采用粒子群算法求解多相平衡下各相中组分的摩尔数；所述粒子群算法包括：

第一步：设置粒子群优化算法的参数和要求的精度，前述粒子群优化算法的参数包括粒子数，迭代次数上限，以及速度由以邻域最优值过渡到全局最优值的迭代分界点；同时获取约束矩阵A以及b；

第二步：获取种群中粒子位置区间范围；

第三步：在粒子位置区间范围内随机初始化粒子群的位置信息，使得粒子的位置位于Ax＝b的超平面内，同时把两次随机初始化的位置信息相减赋给对应种群粒子的速度；

第四步：根据迭代的更新机制，在所述迭代分界点之前利用邻域最优值对速度进行更新，在所述迭代分界点及之后利用全局最优值对速度进行更行；以及依据当前的速度获取更新步长用以更新粒子的位置信息，并更新种群及粒子的历史最优位置；

第五步：检验迭代次数是否达到最初给定值，如达到，种群更新结束，否则返回第四步继续运行，得到平衡时各相中各组分的摩尔数。

本发明具有以下有益效果：

将粒子群算法用于冶金过程多相平衡组分预测，种群中粒子速度的每一维分量是同步进行更新的，这种机制迭代的粒子能够一直保持在可行域中进行解的探索，可以高效实现优化问题中全局最优解的求解。而且，粒子速度的更新是分段进行的；在迭代分界点之前，例如，前90％迭代周期内，利用邻域最优值对速度进行更新，增大了种群中粒子的多样性；而在迭代分界点及之后，例如，后10％迭代周期内，利用全局最优值对速度进行更行，可以使得种群中的粒子尽快的收敛到全局最优解附近。藉此，可使得种群在迭代初期不易早熟，而在迭代末期能够快速收敛至全局最优解。

下面将参照附图，对本发明作进一步详细的说明。

附图说明

构成本申请的一部分的附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：

图1是本发明实施例1公开的复杂冶金过程模拟计算方法流程图；

图2是本发明实施例1公开的粒子群算法的流程图；

图3是本发明实施例1公开的具体运算例中矩阵A的部分示意图；

图4是本发明实施例1公开的具体运算例中b的部分示意图；

图5是本发明实施例1公开的具体运算例中的邻域拓扑结构图；其中图5(a)种群环形拓扑结构示意图，图5(b)为邻域大小为2种群拓扑结构示意图；

图6是本发明实施例2公开的迭代次数与适应度值的关系图；

图7是本发明实施例3公开的迭代次数与适应度值的关系图。

具体实施方式

以下结合附图对本发明的实施例进行详细说明，但是本发明可以由权利要求限定和覆盖的多种不同方式实施。

实施例1

本实施例以铜冶炼为例，公开一种复杂冶金过程模拟计算方法，如图1所示，包括：

步骤S1、以反应体系总的吉布斯自由能函数为该数学模型的目标函数，以输入和输出冶炼过程体系中各种元素的质量相等为约束条件，建立多相平衡数学模型。

在铜冶炼工艺过程中，其反应体系的总吉布斯自由能可由下式(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分别为反应体系中总吉布斯自由能，以及铜锍相，炉渣相和气相吉布斯自由能。而式(2)中，N_p和N_c分别表示体系总相数和各相中组分数，n_ij为j相中组分i的摩尔量，μ_ij表示在反应体系条件下j相中组分i的化学势。ΔG_ij ^o表示在系统温度、标准压强下j相中组分i的吉布斯生成自由能，R是理想气体常数，T表示开尔文温度，f_ij是在j相中组分i的偏逸度，

是j相中组分i在参考状态下的逸度。

受质量守恒和组分摩尔量为非负的约束，这种约束条件在冶金化工领域常见。在铜火法冶炼工艺过程中，由于CaO、MgO、Al₂O₃、SiO₂这类惰性组分直接全部进入渣相，而不参与反应，因而炉渣中这部分组分的摩尔量可以由其进入量而直接得到。而混合矿中所带入的水分，全部转化为水蒸气而进入气相中，因而气相中的水蒸气摩尔量也可以直接得到。同时气相中的氮气呈惰性物质，可以从进入反应体系中的富氧计算得到。其他的元素如Cu、Fe、S、O、As、Sb、Bi、Pb、Zn等由于平衡时候在三相中均有对应的组分存在，因而其质量守恒条件满足，进入体系的元素总量等于其在三相中分配元素的总量。其质量守恒表达式可以用下式(3)表示：

Ax＝b (3)

其中A为由对应组分化学式组成的原子矩阵；x为达到平衡时，由三相中各组分的摩尔量而组成的列向量，即为最终需要求解的变量；而b为由进入体系中元素的总量而组成的列向量。同时由于x中各组分的摩尔量不允许出现负数，所以必须同时保证x>0。

藉此，本实施例所建立多相平衡数学模型为：

min f(x)

st.A·x＝b

x＞0。

步骤S2、采用粒子群算法求解多相平衡下各相中组分的摩尔数。

如图2所示，该步骤可细分为：

步骤S21、设置粒子群优化算法的参数和要求的精度，前述粒子群优化算法的参数包括粒子数，迭代次数上限，以及速度由以邻域最优值过渡到全局最优值的迭代分界点；同时获取约束矩阵A以及b。

步骤S22、获取种群中粒子位置区间范围。换言之，即确定变量x中第j维分量的最大上限值，可以通过下式确定：

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

其中，b_i为体系中相应元素的含量，A_ij为矩阵A中相应的原子系数。

步骤S23、在粒子位置区间范围内随机初始化粒子群的位置信息，使得粒子的位置位于Ax＝b的超平面内，同时把两次随机初始化的位置信息相减赋给对应种群粒子的速度。

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

从矩阵A的n列中选出线性无关的m列，用B表示这m阶方阵，用C表示A中剩下的(n-m)列而组成的m×(n-m)子矩阵，对应的变量x可以相应的分解为x＝[x_B；x_C]，则可以转化为：

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于区间[0,x_Cmax]中，对于求解出的基变量x_B，检验其中x_B的每一维分量是否均大于0，如果是，则x＝[x_B；x_C]即成功初始化在约束超平面Ax＝b内；如果x_B中每一维分量并不都大于0，则需要重新随机赋值x_C于区间[0,x_Cmax]中直至求解出的x_B中每一维分量都大于0；其中，x_Cmax为粒子变量各相应维度的最大上限值。

步骤S24、根据迭代的更新机制，在迭代分界点之前利用邻域最优值对速度进行更新，在迭代分界点及之后利用全局最优值对速度进行更行；以及依据当前的速度获取更新步长用以更新粒子的位置信息，并更新种群及粒子的历史最优位置。

优选的，在迭代更新时，还可以对种群中的粒子增加位置和速度扰动。例如：

位置扰动是通过三个粒子位置向量线性组合生成一个临时的粒子位置pbest_temp实现的，计算公式如下：

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

其中两个粒子位置向量pbest_rand1和pbest_rand2是从当前粒子最优位置池中随机选择，r是区间(0,1)之间的随机数，扰动位置pbest_temp会同相应的当前位置pbest_i进行比较，如果扰动位置pbest_temp相比当前位置pbest_i在种群中有更好的适应度值，则粒子当前的最优位置pbest_i就会被pbest_temp替换。

速度扰动是当前速度通过一个称为速度状态矩阵进行线性变化生成一个临时的粒子速度v_temp实现的，如下式所示：

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

速度状态矩阵v中的行是由从当前粒子速度池中随机选取出的m个粒子速度向量组成；在临时速度基础上，更新当前粒子的位置产生一个新的pbest_temp如下式：

pbest_temp＝pbest_i+λ·v_temp (9)

最后，比较当前粒子的最优位置与新产生的pbest_temp在种群中的适应度值，如果pbest_temp较优，则用pbest_temp替代原先的最优位置pbest_i；其中λ为步长。

其中，在该步骤S24所进行种群中的粒子i更新位置的最大步长λ_{i max}，通过下式获取：

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

其中ε为计算精度，值为用来避免0作为除数但约等于0的正数，x_jmax为变量x中第j维分量的最大上限值，x_ij为粒子i的j维分量的位置，v_ij为粒子i的j维分量的速度。

换言之：上式(9)中，i表示第i个粒子，j表示该粒子的第j维。ε为计算精度用来避免0作为除数，其值约等于0的正数，其中：

当v>ε(其含义可视为v>0)时，视为粒子向正方向运动，其边界为x_i(t+1)＝x_max；

当v<-ε(其含义可视为v<0)时，视为粒子向负方向运动，其边界为x_i(t+1)＝0；

当ε<v<-ε(其含义可视为v＝0)时，视为粒子不运动，此时不需对粒子步长进行约束。

优选的，本实施例中，种群中粒子i的步长系数λ_i按如下公式进行更新：

在上述迭代更新步骤中，所谓“适应度”是将求解出的粒子代入反应体系总吉布斯自由f(x)中，评估标准为：f(x)越小，其适应度越好。

步骤S25、检验迭代次数是否达到最初给定值，如达到，种群更新结束，否则返回步骤S24继续运行，得到平衡时各相中各组分的摩尔数。

为便于本领域技术人员对上述过程做进一步理解，本发明以具体的运算例对上述过程进行复述：

在铜冶炼工程中，x为Cu₂S、Cu、FeS、FeO、Fe₃O₄、Pb、PbS、ZnS、As、Sb、Bi，FeO、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等这些物质在熔炼过程中的含量。其中，铜锍相中主要包括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等组分。

系数矩阵A为上述物质组成的系数矩阵，截图的部分矩阵A如图3所示，其中，矩阵的行是某一元素在上述物质中的原子系数(即脚标)，以金属Cu为例，它在A中表示为：[2 12]。

如图4所示，b为体系中某元素的含量(物质的量)，可根据投入物料计算得到；以金属Cu为例，在加料量66吨，铜含量24.352％时，体系中Cu：66*1000*1000/64＝253107.40157480314mol。

由于f(x)函数对应的等式方程为公式(2)，而公式(2)及公式(3)组成的多相平衡数学模型是欠定方程组，理论上有很多解，这里就是利用粒子群算法求解该方程，使其解满足方程(3)及x>0的同时使公式(2)达到最小值。粒子的位置即为方程的解，也即为体系中物质的含量。

由于所有元素组成的Ax＝b太复杂，这里选取元素铜、铁为例说明：

体系中含铜的物质有三种：Cu₂S、Cu、Cu₂O，设其含量为x₁、x₂、x₃，(其他含铜化合物含量少，忽略不计)

体系中含铁的物质有三种：FeS、FeO、Fe₃O₄,设其含量为x₄、x₅、x₆,(其他含铁化合物含量少，忽略不计)

其Ax＝b可以表示为：

将上述方程(12)展开，得：

对于公式(a)，要想使所有的x均大于0，则每个x必有一个最大值，该最大值必须小于常数项与系数的比值，以x₁为例，其最大值为253107.40/2，因为如果x₁大于该值，其他x必须取负值才能满足上述等式方程。

上述是比较特殊的情况，矩阵A中每一列只有一个值，若矩阵A中每一列有两个或更多的值，如下：

将上述方程展开，得

在该方程组内，如按上述规则取x₁的最大值，会对应两个值(253107.40/2,315385.71/5)，则x₁取最小的一个(315385.71/5)作为最大值，才能使所有的x取值非负。

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

上述式(15)及式(16)可用于解释上述步骤S22中的获取种群中粒子位置区间范围。下面对粒子群的初始化做进一步描述：

为保证计算的速度和精度，在计算过程中会使用若干粒子(数量与实际问题的复杂程度有关)组成的粒子群，由于每个粒子的初始化和迭代过程类似，这里以一个粒子(假设问题只考虑Fe、Cu)为例，该粒子共有6维，每一维即为一个物种(x_i)。

从上式(13)中的矩阵A中取出m(m＝2)线性无关列，构成基向量

剩余的构成非基向量

从新调整

矩阵

使其仍满足

由线性代数知识，在上述范围内任意给定x_C的值，可以求得:

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

其中，x_Cmax为粒子变量各相应维度的最大上限值。通过上式(17)中随机赋值非基变量x_C于区间[0,x_Cmax]中，对于求解出的基变量x_B，检验其中x_B的每一维分量是否均大于0，如果是，则x＝[x_B；x_C]即成功初始化在约束超平面Ax＝b内；如果x_B中每一维分量并不都大于0，则需要重新随机赋值x_C于区间[0,x_Cmax]中直至求解出的x_B中每一维分量都大于0；即得到了问题的初始解(一个粒子的初始位置)。粒子在运动过程中，两次位置之差即为粒子运动的速度v，因此，可以对每个粒子初始化两次，将两次位置之差记为粒子的初始化速度。

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

由于上述初始化满足Ax＝b，即为方程组的解，但是并不一定能满足目标函数值达到最小，因此需要对粒子群的位置进行更新，使其满足Ax＝b，且x>0的前提下，使目标函数值f(x)最小。为此，关于上述步骤S24中的例子速度及位置更新以及扰动等，通过下述(一)至(四)详述如下：

(一)、种群迭代步长求解

为了避免种群粒子在下次速度更新完之后越过可行解搜索域，需要对初始速度进行相应的速度大小限制，通过增加步长系数λ即可满足粒子每次速度更新后不会越过可行解搜索域。种群中粒子i的最大步长系数λ_{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}，相应的实验表明种群的局部搜索性能会逐渐减小，因而在本实施例，当粒子的λ_{i max}>1时候，这部分粒子的更新步长系数重新被置为1，而其他粒子保持不变。因而种群中粒子i的步长系数λ_i按如下公式(11)进行更新。

(二)、种群拓扑结构设计

实际过程中，对于复杂难优化的问题，要求种群有好的全局搜索能力，而当种群定位到最优解附近位置时，要求加强局部探索能力，即获取精度更高的最优解。好的种群拓扑结构能够有效的保证搜索初期在全局范围内避免种群早熟的探寻能力，以及搜索后期在局部范围内更精确的探寻能力。

为此，如附图5(a)所示，种群拓扑结构可采用环形结构；本实施例设计了一种名为单链环形邻域拓扑结构的粒子群，一种典型的邻域大小为2(邻域大小为所考虑的周围粒子的数量，不代表距离)的种群拓扑结构如附图5(b)所示，对于种群中粒子k，邻域由种群中粒子编号为k+1，k-2，k+3，k-4…k+n-1和k-n的N个粒子组成；这种拓扑结构，能够有效避免编号相邻粒子间的相互吸引。通过实施上述选择邻域的策略，种群中粒子i的邻域最佳粒子记为lbest_i。

(三)、种群全局搜索与局部探索能力平衡机制

从全局的视角出发平衡种群的全局搜索能力和局部探索能力，对种群的权重因子影响和速度更新策略进行了充分的考虑。一般，迭代初期需要保持种群有较强的全局搜索能力，种群的粒子能够尽可能分散在搜索域的各个角落；而迭代后期需要保持种群有较强的局部探索能力，种群的粒子能够尽快收敛到全局最优值附近，满足探索更高精度的优化解。

粒子群算法中，高的权重因子就意味着粒子能够快速的在搜索空间内移动，保证粒子能够达到搜索空间的任意位置。因而权重因子在迭代初期应该保持高的水平，而随着迭代次数的增加，权重因子w逐渐降低到零。在本发明中的粒子群优化算法中，权重因子的更新遵循下面的公式(19)：

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

在上述公式中，MaxIter是种群的总迭代次数，k是当前的迭代次数。藉此，权重因子随着迭代次数的变化呈现一条S曲线进行变化。

在本实施例中，粒子速度的更新是分段进行的。例如：在前90％迭代周期内，种群中粒子的速度更新遵循公式(20)，为了增大种群中粒子的多样性；而在后10％迭代周期内，种群中粒子的速度更新遵循公式(21)，为了使得种群中的粒子尽快的收敛到最优解附近。因此上述的策略能够使得种群在迭代初期不易早熟，而在迭代末期能够快速收敛至全局最优解。种群中粒子速度的每一维分量是同步进行更新的，这种机制迭代的粒子能够一直保持在可行域中进行解的探索。

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为粒子当前的最优位置，lbest_i为粒子邻域内的最优位置，gbest_i为粒子群内当前的最优位置，r₁和r₂是区间(0,1)之间两相互独立的随机数。

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

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

(四)、种群粒子位置和速度扰动

对于当前粒子i最优位置pbest_i以及速度v_i添加扰动能够在保持种群多样性的同时，引导粒子飞向更优的位置而不破坏种群自身的组织结构。

种群中粒子位置扰动是通过三个粒子位置向量线性组合生成一个临时的粒子位置pbest_temp实现的，如式(7)所示。其中两个粒子位置向量是从的当前粒子最优位置池中随机选择，r是区间(0,1)之间的随机数。通过对当前位置pbest_i添加扰动，使得产生的pbest_temp仍然满足A·pbest_temp＝b。之后扰动位置pbest_temp会同相应的当前位置pbest_i进行比较，如果扰动位置pbest_temp相比当前位置pbest_i在种群中有更好的适应度值，则粒子当前的最优位置pbest_i就会被pbest_temp替换。

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

种群中对粒子速度的扰动是把当前速度通过一个称为速度状态矩阵进行线性变化生成一个临时的粒子速度v_temp实现的，如式8)所示。

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

矩阵v即为速度状态矩阵，是一个m×m的方阵，m即为粒子群的大小。速度状态矩阵v中的行是由从当前粒子速度池中随机选取出的m个粒子速度向量组成。因而当前的粒子速度v_i在速度状态矩阵v转化下，能够向任何方向进行偏移，种群的局部搜索能力能够得到大大的提高。而经过此扰动，生成的v_temp仍然能够满足A·v_temp＝0条件。在此基础上，根据式(10)和(11)确定每次迭代的步长，更新当前粒子的位置产生一个新的pbest_temp如式(9)所示。最后，比较当前粒子的最优位置与新产生的pbest_temp在种群中的适应度值，如果pbest_temp较优，则用pbest_temp替代原先的最优位置pbest_i。

pbest_temp＝pbest_i+λ·v_temp (9)

综上，本发明适用于求解带超高维数线性约束，且非凸目标函数优化问题的改进粒子群优化算法(HLPSO)的主旨如下：最先需要获取原子矩阵A和列向量b，以及算法的相关设置参数。由于约束条件变量x非负的要求，同时A和b中的成员均为非负，对于变量x来说自然存在上界。保持变量x不超出上界保证了后续的迭代操作均约束在Ax＝b超平面内。之后即为最重要的两个步骤：种群初始化和种群过程迭代。

HLPSO算法中种群的初始化包括粒子位置初始化和速度初始化。由于粒子初始位置是随机生成的，因而有一定的概率初始化失败，直至初始化所有粒子的位置x在0与上界之间的区间。而粒子的速度初始化即为两次粒子随机初始化的位置之差。在所有粒子均初始化成功后，选择种群中所有粒子位置中的最优位置为全局最优位置。上述种群初始化成功后，即进入种群过程迭代。

种群迭代过程即粒子学习过程，逐步靠近全局最优值的过程。每次迭代过程中，获取种群全局及粒子邻域最优位置用于粒子速度更新，而速度更新过程中需要对速度做相应的限制，以免在位置更新时候跃出可行搜索区域。之后更新种群及粒子的历史最优位置，在此基础上从粒子池中随机选择粒子数据添加对当前粒子的位置和速度的扰动，以避免迭代前期陷入局部最优值，同时在迭代后期能够加强粒子局部探索能力。

迭代次数达到后，种群更新结束，而最后种群中的历史全局最优位置极可能为所求解问题的全局最优值。

实施例2：

为了测试算法的性能，从文献中选取了如下式(23)的测试函数，其函数形式同本设计铜冶炼过程多相平衡数学模型类型相类似，具有一定的代表性。

其中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.设计提出的，ISRES算法是由Runarsson和Yao提出的。其与HLPSO针对上述测试函数的测试结果比较如下表1所示。

表1：

文献中给出的精确解为：

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

而使用HLPSO每次都能收敛到上述最优解附近，十次测试的平均最优解为：x＝(0.0407，0.1477，0.7831，0.0014，0.4853，6.959e-04，0.0274，0.018，0.0373，0.0968)，由此可见与文献中给出的精确解几乎一样，说明HLPSO算法的是可行有效的。

由上述比较可知，HLPSO算法在处理约束方面的内在特性，种群能够始终保持在可行域内，因而搜索效率相比于COPSO和ISRES性能大大提高，其稳定性也好于其他两种算法，由图6可以知道，HLPSO算法很快的就能够搜索到全局最优值附近，在总计100次迭代内，大约40次迭代后就收敛到了全局最优值，收敛速度块，具有较高的搜索效率。因而HLPSO算法是适合针对本发明的所提出的铜冶炼过程多相平衡数学模型的求解。

实施例3：

在国内某铜冶炼厂，采集了相关生产实践过程中的基础参数，建立多相平衡数学模型。

通过生产调研，冰铜和炉渣相的温度约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等组分。

在上述基础上通过热力学多相平衡理论建立数学模型，确定了优化目标函数以及约束条件，并对其进行计算求取最终平衡时候各相中的组分含量。

从国内某以铜冶炼厂获取了其在2014年稳定工况下某段时间内的入炉混合铜精矿成分表及其对应的操作工艺条件。其中表2为入炉铜精矿的成分表，而工艺操作过程参数如表3所示。

表2：

续

表3：

在上述生产工况下，选取处理得到的工业生产数据同模拟计算的结果进行了对比，如下表4和表5所示。

表4：

表5：

通过以上模拟计算数据和工业生产数据对比分析结果可知，所建立的铜冶炼过程多相平衡模型是可靠的，可用于仿真预测铜冶炼工艺过程中多元素的分配行为以及各组分之间的相互关系。

设置HLPSO算法中种群大小为200，而迭代次数为1000，重复进行相同的测试10次。10次重复测试实验的结果如下表6所示。

表6：

由上述测试可知，针对铜冶炼过程多相平衡数学模型，HLPSO算法能够胜任上述模型的计算，计算成功率为90％，几乎每次都能收敛至全局最优值附近。图7说明了HLPSO算法很快就能够搜索到全局最优值附近，在总计1000次迭代内，大约200次迭代后就收敛到了全局最优值，收敛速度块，具有较高的搜索效率。因而HLPSO算法能够适用于本发明所提出的铜冶炼过程多相平衡数学模型的求解，且HLPSO算法是可靠的。

实施例4

与上述实施例相对应的，本发明实施例还公开一种复杂冶金过程模拟计算系统，至少包括：

多相平衡建模模块，用于以反应体系总的吉布斯自由能函数为该数学模型的目标函数，以输入和输出冶炼过程体系中各种元素的质量相等为约束条件，建立多相平衡数学模型。

如上述实施例1，该多相平衡建模模块所建立的多相平衡数学模型具体为：

min f(x)

st.A·x＝b

x＞0

上述模型中，矩阵A是由各相中各组分系数组成的原子系数矩阵，对应的b是由进入反应体系中各元素的总摩尔量组成的列向量，x为各相中各组分的摩尔数，f(x)为反应体系总吉布斯自由能；

算法处理模块，用于采用粒子群算法求解多相平衡下各相中组分的摩尔数。所述粒子群算法包括：

第一步：设置粒子群优化算法的参数和要求的精度，前述粒子群优化算法的参数包括粒子数，迭代次数上限，以及速度由以邻域最优值过渡到全局最优值的迭代分界点；同时获取约束矩阵A以及b；

第二步：获取种群中粒子位置区间范围；

第三步：在粒子位置区间范围内随机初始化粒子群的位置信息，使得粒子的位置位于Ax＝b的超平面内，同时把两次随机初始化的位置信息相减赋给对应种群粒子的速度；

第四步：根据迭代的更新机制，在所述迭代分界点之前利用邻域最优值对速度进行更新，在所述迭代分界点及之后利用全局最优值对速度进行更行；以及依据当前的速度获取更新步长用以更新粒子的位置信息，并更新种群及粒子的历史最优位置；

第五步：检验迭代次数是否达到最初给定值，如达到，种群更新结束，否则返回第四步继续运行，得到平衡时各相中各组分的摩尔数。

其中，本实施例中算法处理模块的具体计算请参照上述实施例1，在此不做赘述。综上，本实施例所公开的复杂冶金过程模拟计算方法及系统，可以高效实现冶金过程多相平衡组分预测等优化问题中全局最优解的求解。具体的：将粒子群算法用于冶金过程多相平衡组分预测，种群中粒子速度的每一维分量是同步进行更新的，这种机制迭代的粒子能够一直保持在可行域中进行解的探索，可以高效实现优化问题中全局最优解的求解。而且，粒子速度的更新是分段进行的；在迭代分界点之前，例如，前90％迭代周期内，利用邻域最优值对速度进行更新，增大了种群中粒子的多样性；而在迭代分界点及之后，例如，后10％迭代周期内，利用全局最优值对速度进行更行，可以使得种群中的粒子尽快的收敛到全局最优解附近。藉此，可使得种群在迭代初期不易早熟，而在迭代末期能够快速收敛至全局最优解。

此外，优选的，当本实施例的上述方法及系统应用于冶炼铜时，还可以进一步将多相平衡数学模型与机械夹杂方程结合起来，对主要成分进行机械夹杂的修正；其中：

铜锍在炉渣中的夹杂率计算公式为：

炉渣在铜锍相中的夹杂率计算公式为：

铜锍相中夹杂炉渣量的计算公式为：

炉渣相中夹杂铜锍量的计算公式为：

其中，M_slag和M_matte分别表示理论计算平衡时炉渣相和铜锍相的总质量，M_slag和M_matte分别是进入铜锍相中的炉渣质量和进入炉渣相中的铜锍质量，

和

分别是炉渣在铜锍相中的夹杂系数以及铜锍在炉渣相中的夹杂系数。

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种复杂冶金过程模拟计算方法，其特征在于，包括：

获取稳定工况下设定时间段内的入炉矿物成分以及工艺操作参数；

以在所述入炉矿物成分以及工艺操作参数条件下的反应体系总的吉布斯自由能函数为数学模型的目标函数，以输入和输出冶炼过程体系中各种元素的质量相等为约束条件，建立多相平衡数学模型：

min f(x)

st.A·x＝b

x>0

矩阵A是由各相中各组分系数组成的原子系数矩阵，对应的b是由进入反应体系中各元素的总摩尔量组成的列向量，x为各相中各组分的摩尔数，f(x)为反应体系总吉布斯自由能；

采用粒子群算法求解多相平衡下各相中组分的摩尔数；所述粒子群算法包括：

第一步：设置粒子群算法的参数和要求的精度，前述粒子群算法的参数包括粒子数，迭代次数上限，以及速度由以邻域最优值过渡到全局最优值的迭代分界点；同时获取矩阵A以及b；

第二步：获取种群中粒子位置区间范围；

第三步：在粒子位置区间范围内随机初始化粒子群的位置信息，使得粒子的位置位于A·x＝b的超平面内，同时把两次随机初始化的位置信息相减赋给对应种群粒子的速度；

第四步：根据迭代的更新机制，在所述迭代分界点之前利用邻域最优值对速度进行更新，在所述迭代分界点及之后利用全局最优值对速度进行更行；以及依据当前的速度获取更新步长用以更新粒子的位置信息，并更新种群及粒子的历史最优位置；

第五步：检验迭代次数是否达到最初给定值，如达到，种群更新结束，否则返回第四步继续运行，得到平衡时各相中各组分的摩尔数；

预测得到铜冶炼工艺过程中多元素的分配行为以及各组分之间的相互关系。

2.根据权利要求1所述的复杂冶金过程模拟计算方法，其特征在于，所述获取种群中粒子位置区间范围即确定变量x中第j维分量的最大上限值，通过下式确定：

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

其中，b_i为体系中相应元素的含量，A_ij为矩阵A中相应的原子系数。

3.根据权利要求1所述的复杂冶金过程模拟计算方法，其特征在于，初始化位置和速度信息包括：

从矩阵A的n列中选出线性无关的m列，用B表示这m阶方阵，用C表示A中剩下的(n-m)列而组成的m×(n-m)子矩阵，对应的变量x可以相应的分解为x＝[x_B；x_C]，则转化为：

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

则x_B即通过下式得到：

x_B＝B^-1·b-B^-1·C·x_C

通过随机赋值非基变量x_C于区间[0,x_Cmax]中，对于求解出的基变量x_B，检验其中x_B的每一维分量是否均大于0，如果是，则x＝[x_B；x_C]即成功初始化在约束超平面A·x＝b内；如果x_B中每一维分量并不都大于0，则需要重新随机赋值x_C于区间[0,x_Cmax]中直至求解出的x_B中每一维分量都大于0；

其中，x_Cmax为粒子变量各相应维度的最大上限值。

4.根据权利要求1所述的复杂冶金过程模拟计算方法，其特征在于，种群中的粒子i更新位置的最大步长λ_imax，通过下式获取：

λ_imax＝min(x_jmax-x_ij/v_ij|v_ij>ε,-x_ij/v_ij|v_ij<-ε,0|-ε<v_ij<ε)；

其中ε为计算精度，值为用来避免0作为除数但约等于0的正数，x_jmax为变量x中第j维分量的最大上限值，x_ij为粒子i的j维分量的位置，v_ij为粒子i的j维分量的速度；以及种群中粒子i的步长系数λ_i按如下公式进行更新：

5.根据权利要求1所述的复杂冶金过程模拟计算方法，其特征在于，在迭代更新速度和位置信息时，还包括：

对种群中的粒子增加位置和速度扰动；和/或

采用邻域大小大于或等于2的单链环形邻域拓扑结构。

6.根据权利要求5所述的复杂冶金过程模拟计算方法，其特征在于，位置扰动是通过三个粒子位置向量线性组合生成一个临时的粒子位置pbest_temp实现的，计算公式如下：

pbest_temp＝pbest_i+r·(pbset_rand1-pbset_rand2)

其中两个粒子位置向量pbest_rand1和pbest_rand2是从当前粒子最优位置池中随机选择，r是区间(0,1)之间的随机数，粒子位置pbest_temp会同相应的当前位置pbest_i进行比较，如果粒子位置pbest_temp相比当前位置pbest_i在种群中有更好的适应度值，则粒子的当前位置pbest_i就会被pbest_temp替换。

7.根据权利要求6所述的复杂冶金过程模拟计算方法，其特征在于，速度扰动是当前速度通过一个称为速度状态矩阵进行线性变化生成一个临时的粒子速度v_temp实现的，如下式所示：

v_temp＝v_i·v/||v||

速度状态矩阵v中的行是由从当前粒子速度池中随机选取出的m个粒子速度向量组成；在临时速度基础上，更新当前粒子的位置产生一个新的pbest_temp如下式：

pbest_temp＝pbest_i+λ·v_temp

最后，比较当前粒子的最优位置与新产生的pbest_temp在种群中的适应度值，如果pbest_temp较优，则用pbest_temp替代原先的最优位置pbest_i；其中λ为步长。

8.根据权利要求6或7所述的复杂冶金过程模拟计算方法，其特征在于，所述适应度值是将求解出的粒子代入反应体系总吉布斯自由能f(x)中，评估标准为：

f(x)越小，其适应度值越好。

9.根据权利要求1所述的复杂冶金过程模拟计算方法，其特征在于，当将该方法应用于冶炼铜时，还包括：

将多相平衡数学模型与机械夹杂方程结合起来，对各组分进行机械夹杂的修正；其中：

铜锍在炉渣中的夹杂率计算公式为：

炉渣在铜锍相中的夹杂率计算公式为：

铜锍相中夹杂炉渣量的计算公式为：

炉渣相中夹杂铜锍量的计算公式为：

其中，M_slag和M_matte分别表示理论计算平衡时炉渣相和铜锍相的总质量，

和

分别是进入铜锍相中的炉渣质量和进入炉渣相中的铜锍质量，

和

分别是炉渣在铜锍相中的夹杂系数以及铜锍在炉渣相中的夹杂系数。

10.一种用于执行权利要求1至9任一所述复杂冶金过程模拟计算方法的系统，其特征在于，至少包括：

输入模块，用于获取稳定工况下设定时间段内的入炉矿物成分以及工艺操作参数；

多相平衡建模模块，用于以在所述入炉矿物成分以及工艺操作参数条件下的反应体系总的吉布斯自由能函数为数学模型的目标函数，以输入和输出冶炼过程体系中各种元素的质量相等为约束条件，建立多相平衡数学模型：

min f(x)

st.A·x＝b

x>0

矩阵A是由各相中各组分系数组成的原子系数矩阵，对应的b是由进入反应体系中各元素的总摩尔量组成的列向量，x为各相中各组分的摩尔数，f(x)为反应体系总吉布斯自由能；

算法处理模块，用于采用粒子群算法求解多相平衡下各相中组分的摩尔数；所述粒子群算法包括：

第一步：设置粒子群算法的参数和要求的精度，前述粒子群算法的参数包括粒子数，迭代次数上限，以及速度由以邻域最优值过渡到全局最优值的迭代分界点；同时获取矩阵A以及b；

第二步：获取种群中粒子位置区间范围；

第三步：在粒子位置区间范围内随机初始化粒子群的位置信息，使得粒子的位置位于A·x＝b的超平面内，同时把两次随机初始化的位置信息相减赋给对应种群粒子的速度；

第四步：根据迭代的更新机制，在所述迭代分界点之前利用邻域最优值对速度进行更新，在所述迭代分界点及之后利用全局最优值对速度进行更行；以及依据当前的速度获取更新步长用以更新粒子的位置信息，并更新种群及粒子的历史最优位置；

第五步：检验迭代次数是否达到最初给定值，如达到，种群更新结束，否则返回第四步继续运行，得到平衡时各相中各组分的摩尔数；

预测得到铜冶炼工艺过程中多元素的分配行为以及各组分之间的相互关系。

Images