CN102982376B - 一种基于遗传计算的二维泊松方程求解优化方法 - Google Patents

Abstract

一种基于遗传计算的二维泊松方程快速求解方法，包括以下步骤：（1）采用遗传计算对松弛因子进行全局寻优，适应度函数建模；（2）初始化种群，对种群进行优胜劣汰的筛选；（3）新个体由父个体的线性插值及非均匀变异产生，交叉概率和变异概率根据自适应遗传算法进行计算；（4）判断收敛性，算法收敛时适应度最大值对应的个体即为最佳松弛因子；（5）若算法收敛，则进行并行超松弛迭代计算，实现二维泊松方程的快速计算；若算法不收敛，则返回步骤（2），继续通过截断选择法与稳态繁殖法相结合对种群进行优胜劣汰的筛选。本发明迭代次数少、求解精度高、计算速度快。

Description

一种基于遗传计算的二维泊松方程求解优化方法

技术领域

本发明涉及信号处理领域，更具体的说是涉及一种基于遗传计算的二维泊松方程求解优化方法。

背景技术

泊松方程的常用解法有格林函数法、分离变量法、有限差分法、迭代算法等。格林函数可以将微分方程边值问题转化为积分方程问题，但对于有限域的泊松方程，因为找不到其对应的格林函数，故该法求解比较困难。分离变量法是求解数学物理方程中应用最广泛的一种方法，但该法在应用时坐标系的选择有一定限制，当所求解模型的边界面与某坐标系统的坐标面相吻合，或者至少分段地与坐标面相吻合时才能使用。有限差分法是最早且很好的求解方法，对于研究抛物型和椭圆型问题，受到人们的关注和重视。但有限差分方法对于椭圆型问题的逼近往往需要求解较大的稀疏矩阵，数据处理也很复杂。并行超松弛迭代算法因其有明显的并行性，能大大提高计算效率、节省计算时间、减少迭代次数，由于采用了并行技术，计算机各处理器间通信与计算时间重叠，获得了较为理想的加速效率，缺点是最佳松弛因子选择困难。在一般的情况下，并行超松弛迭代算法最佳收敛因子只能凭借经验取值，因此如何快速选取最佳因子成为并行超松弛迭代算法的关键。

发明内容

有鉴于此，为了解决并行超松弛迭代算法松弛因子选择困难而影响计算速度问题，本发明提供一种能有效减少迭代次数，提高算法效率，节省计算时间的基于遗传计算的二维泊松方程求解优化方法，将遗传算法与并行超松弛迭代算法相结合，以加快二维泊松方程的求解速度以及提高计算精度。

本发明的技术方案包括以下步骤：

(1)采用现有遗传算法技术对松弛因子进行全局寻优，适应度函数建模为与迭代次数N以及收敛精度有关的多目标适应度函数，式中，u_i,j表示电场、磁场或温度场中点i,j处的位函数，K表示当前迭代数；

(2)初始化种群，采用现有截断选择法与稳态繁殖法相结合对种群进行优胜劣汰的筛选，只保留精英个体，提高种群的多样性；

(3)新个体由父个体的线性插值及非均匀变异产生，交叉概率和变异概率根据自适应遗传算法进行计算；

(4)判断收敛性，遗传算法的收敛条件是迭代次数超过300或者最大适应度连续3代变化都小于10^-10，算法收敛时适应度最大值对应的个体即为最佳松弛因子；

(5)若算法收敛，则选用五台处理能力相同的PC机作为硬件平台，一台作为主机，其余四台作为从机，主机与从机通信，从机之间互不干扰；将遗传计算得到的最佳松弛因子作为并行超松弛迭代算法的松弛因子，进行并行超松弛迭代计算，实现二维泊松方程的快速计算。

若算法不收敛，则返回步骤(2)，继续通过截断选择法与稳态繁殖法相结合对种群进行优胜劣汰的筛选。

所述步骤(1)中，为初始化种群，设置种群大小M＝20，染色体长度为15，进化代数H＝10，变异算子b＝2，截断阈值T＝10％。多目标适应度函数为：

$J_{fit}=\frac{1}{1+c_{1}N+c_{2}max|u_{i,j}^{\left(N\right)}\left(\mathrm{\ω}\right)-u_{i,j}^{(N-1)}\left(\mathrm{\ω}\right)|}---\left(1a\right)$

$\frac{c_1}{c_2}=\frac{{\mathrm{maxJ}}_n^2}{{\mathrm{maxN}}_n},(n=1,2,...,M)---\left(2a\right)$

式中c1、c2为正加权系数，满足c1+c2＝1，N为迭代次数，u_i,j表示电场、磁场或温度场中点i,j处的位函数，ω表示松弛因子，表示当前代的最大误差；maxN_n表示当前代迭代次数的最大值，M为正整数。

所述步骤(3)中线性插值方法的计算公式为：

$\left\{\begin{array}{c}{\mathrm{\ω}}_1^{*}={\mathrm{a\ω}}_1+\left(1-a\right){\mathrm{\ω}}_2\\ {\mathrm{\ω}}_2^{*}=\left(1-a\right){\mathrm{\ω}}_1+{\mathrm{a\ω}}_2\end{array}\right.---\left(3a\right)$

其中a是(0,1)区间内的随机数，ω₁、ω₂为旧个体，为新个体。

所述步骤(3)中非均匀实值变异算子的计算公式为：

${\mathrm{\ω}}^{\′}=\left\{\begin{array}{cc}\mathrm{\ω}+\mathrm{\Δ}\left(t,{\mathrm{\ω}}_{\max }-\mathrm{\ω}\right),& d>0.5\\ \mathrm{\ω}-\mathrm{\Δ}\left(t,\mathrm{\ω}-{\mathrm{\ω}}_{\min }\right),& d\≤0.5\end{array}\right.---\left(4a\right)$

其中ω_max和ω_min分别为ω的上界和下界值，t为当前进化代数，ω为变异前的值，ω'为变异后的值，d为[0,1]区间的随机数，△(t,y)＝y□r□(1-t/W)^b，这里r为[0,1]区间的随机数，W为最大进化代数；b为确定非均匀度的参数；y是中间参数。

所述步骤(3)中交叉概率P_c和变异概率P_m计算公式分别为

$P_c=\left\{\begin{array}{cc}{P}_{c1}-\frac{(P_{c1}-P_{c2})(J^{\&prime;}-J_{avg})}{J_{max}-J_{avg}},& J^{\&prime;}\&GreaterEqual;J_{avg}\\ P_{c1},& J^{\&prime;}<J_{avg}\end{array}\right.---\left(5a\right)$

$P_m=\left\{\begin{array}{cc}{P}_{m1}-\frac{(P_{m1}-P_{m2})(J_{max}-J)}{J_{max}-J_{avg}},& J\&GreaterEqual;J_{avg}\\ P_{m1},& J<J_{avg}\end{array}\right.---\left(6a\right)$

式中J_max为每代群体最大适应度值；J_avg为每代群体平均适应度值；J'为两个交叉个体较大的适应度值；J为变异个体的适应度值；P_c1＝0.9，P_c2＝0.6，P_m1＝0.1，P_m2＝0.001。

一台主机与四台从机构成机群系统，算法通过机群系统实施，主机负责与四台从机进行通信并更新分界线上点的值，四台从机进行分块场域的超松弛迭代计算。

本发明结构简单，迭代次数少，求解精度高；其采用遗传计算对松弛因子进行全局寻优，将遗传计算得到的最佳松弛因子作为并行超松弛迭代算法的松弛因子，进行迭代计算，实现二维泊松方程的快速计算。优点在于构造的多目标适应度函数是一个凸函数，对松弛因子具有全局搜索能力。在保证计算精度的前提下，提高了方程的收敛速度，提高了计算效率，计算速度快。

附图说明

图1为本发明实施例流程框图；

图2为并行超松弛迭代算法格式图；

图3(a)为种群适应度曲线；

图3(b)为图3(a)的局部放大曲线；

图3(c)为图3(b)的局部放大曲线；

图3(d)为每代最大适应度变化曲线；

图3(e)为最佳因子变化曲线；

图3(f)为每代最大误差变化曲线。

具体实施方式

下面结合附图和实施例对本发明作进一步详细的说明。

参照图1，本发明所提出的算法包括以下计算步骤：

(1)首先执行步骤101，种群大小M＝20，染色体长度为15，进化代数H＝10，变异算子b＝2，截断阈值T＝10％，染色体编码方式为实值编码，完成种群的初始化。

(2)执行步骤102，采用截断选择策略，利用多目标适应度函数区分个体的优劣，其表达式为：

$J_{fit}=\frac{1}{1+c_{1}N+c_{2}max|u_{i,j}^{\left(N\right)}\left(\mathrm{\&omega;}\right)-u_{i,j}^{(N-1)}\left(\mathrm{\&omega;}\right)|}---\left(1b\right)$

$\frac{c_1}{c_2}=\frac{{\mathrm{maxJ}}_n^2}{{\mathrm{maxN}}_n},(n=1,2,...,M)---\left(2b\right)$

式中c1、c2为正加权系数，满足c1+c2＝1，N为迭代次数，u_i,j表示电场、磁场或温度场中点i,j处的位函数，ω表示松弛因子，表示当前代的最大误差；maxN_n表示当前代迭代次数的最大值，M为正整数。

采用稳态繁殖法，使当前群体中所有适应度最好的个体尽可能地保留到下一代群体中，其具体操作过程为：首先找出当前群体中所有适应度最高的个体；然后对当前群体进行选择、交叉、变异运算后产生的新群体的个体适应度进行排序；最后用当前群体中所有适应度最高的个体替换新群体中适应度低的个体。按上述方法设计算法的优点主要表现在两个方面：减小了算法的搜索范围，同时交叉和变异运算不会破坏适应度最好的个体；在选择、交叉、变异运算后采用稳态繁殖法可以加快算法的收敛速度，并能保证找到所有的最优解。

(3)执行步骤103，根据交叉概率P_c选择父个体进行线性插值产生新的子个体，其计算公式为：

$\left\{\begin{array}{c}{\mathrm{\&omega;}}_1^{*}={\mathrm{a\&omega;}}_1+\left(1-a\right){\mathrm{\&omega;}}_2\\ {\mathrm{\&omega;}}_2^{*}=\left(1-a\right){\mathrm{\&omega;}}_1+{\mathrm{a\&omega;}}_2\end{array}\right.---\left(3b\right)$

其中a是(0,1)区间内的随机数，ω₁、ω₂为旧个体，为新个体。交叉概率P_c的表达式为：

$P_c=\left\{\begin{array}{cc}{P}_{c1}-\frac{(P_{c1}-P_{c2})(J^{\&prime;}-J_{avg})}{J_{max}-J_{avg}},& J^{\&prime;}\&GreaterEqual;J_{avg}\\ P_{c1},& J^{\&prime;}<J_{avg}\end{array}\right.---\left(4b\right)$

式中J_max为每代群体最大适应度值；J_avg为每代群体平均适应度值；J'为两个交叉个体中较大者对应的适应度值；J为变异个体的适应度值；P_c1＝0.9，P_c2＝0.6，P_m1＝0.1，P_m2＝0.001。

(4)执行步骤104，交叉得到的个体再按变异概率P_m进行非均匀变异，表达式为：

${\mathrm{\&omega;}}^{\&prime;}=\left\{\begin{array}{cc}\mathrm{\&omega;}+\mathrm{\&Delta;}\left(t,{\mathrm{\&omega;}}_{\max }-\mathrm{\&omega;}\right),& d>0.5\\ \mathrm{\&omega;}-\mathrm{\&Delta;}\left(t,\mathrm{\&omega;}-{\mathrm{\&omega;}}_{\min }\right),& d\&le;0.5\end{array}\right.---\left(5b\right)$

其中ω_max和ω_min分别为ω的上界和下界值，t为当前进化代数，ω为变异前的值，ω'为变异后的值，d为[0,1]区间的随机数，△(t,y)＝y□r□(1-t/W)^b，这里r为[0,1]区间的随机数，W为最大进化代数；b为确定非均匀度的参数；y是中间参数。变异概率P_m的表达式为：

$P_m=\left\{\begin{array}{cc}{P}_{m1}-\frac{(P_{m1}-P_{m2})(J_{max}-J)}{J_{max}-J_{avg}},& J\&GreaterEqual;J_{avg}\\ P_{m1},& J<J_{avg}\end{array}\right.---\left(6b\right)$

式中J_max为每代群体最大适应度值；J_avg为每代群体平均适应度值；J'为两个交叉个体较大的适应度值；J为变异个体的适应度值；P_c1＝0.9，P_c2＝0.6，P_m1＝0.1，P_m2＝0.001。

(5)执行步骤105，判断收敛性，遗传算法的收敛条件是迭代次数超过300或者最大适应度连续3代变化都小于10^-10，算法收敛时适应度最大值对应的个体即为最佳松弛因子。如果算法收敛则执行步骤(6)，否则返回执行步骤(2)。

(6)执行步骤106，选用5台处理能力相同的PC机作为硬件平台，一台作为主机，其余四台作为从机，主机与从机通信，从机之间互不干扰。四个并行超松弛迭代计算：

$\left\{\begin{array}{c}{u}_{i,j}^{\left(n+1\right)}=\left(1-\mathrm{\&omega;}\right)u_{i,j}^{\left(n\right)}+\frac{\mathrm{\&omega;}}{4}\left(u_{i+1,j}^{\left(n+1\right)}+u_{i,j+1}^{\left(n+1\right)}+u_{i-1,j}^{\left(n\right)}+u_{i,j-1}^{\left(n\right)}+h^2{f}_{i,j}\right)\\ u_{i,j+1}^{\left(n+1\right)}=\left(1-\mathrm{\&omega;}\right)u_{i,j+1}^{\left(n\right)}+\frac{\mathrm{\&omega;}}{4}\left(u_{i+1,j+1}^{\left(n+1\right)}+u_{i,j+2}^{\left(n\right)}+u_{i-1,j+1}^{\left(n\right)}+u_{i,j}^{\left(n+1\right)}+h^2{f}_{i,j+1}\right)\\ u_{i+1,j}^{\left(n+1\right)}=\left(1-\mathrm{\&omega;}\right)u_{i+1,j}^{\left(n\right)}+\frac{\mathrm{\&omega;}}{4}\left(u_{i+2,j}^{\left(n\right)}+u_{i+1,j+1}^{\left(n+1\right)}+u_{i,j}^{\left(n+1\right)}+u_{i+1,j-1}^{\left(n\right)}+h^2{f}_{i+1,j}\right)\\ u_{i+1,j+1}^{\left(n+1\right)}=\left(1-\mathrm{\&omega;}\right)u_{i+1,j+1}^{\left(n\right)}+\frac{\mathrm{\&omega;}}{4}\left(u_{i+2,j+1}^{\left(n\right)}+u_{i+1,j+2}^{\left(n\right)}+u_{i,j+1}^{\left(n+1\right)}+u_{i+1,j}^{\left(n+1\right)}+h^2{f}_{i+1,j+1}\right)\end{array}\right.---\left(7b\right)$

式中，u_i,j表示电场、磁场或温度场点i,j处的位函数，n表示当前迭代数，ω表示松弛因子，h表示场域的网格划分长度，f表示场域(电场、磁场、温度场等)中二维泊松方程满足的条件。

参照图2，将目标问题所在场域(电场、磁场、温度场等)划分为四个部分，用四台从机对子域A、B、C、D分别进行超松弛迭代，主机负责迭代子区域上的边界点并将新值传送到从机。每次迭代后验证收敛性，即判断泊松方程是否达到精度要求，若达到则停止迭代。

本发明应用实施例：

参照图3，假设某圆筒的温度分布与Z轴无关，并且假设由于感应加热而产生的热源强度为：

$q(r,\mathrm{\&theta;})=\frac{q_0}{2}{r}^2(1+cos2\mathrm{\&theta;})---\left(8b\right)$

式中，q₀为常数，r表示与热源的距离，θ为角度。

此时，圆筒内部的温度u(r,θ)满足二维泊松方程：

${\&dtri;}^{2}u+\frac{q}{k}=0---\left(9b\right)$

式中，q为热源强度，k为常数。

或

$\frac{{\&part;}^{2}u}{\&part;r^2}+\frac{1}{r}\frac{\&part;u}{\&part;r}+\frac{1}{r^2}\frac{{\&part;}^{2}u}{\&part;{\mathrm{\&theta;}}^2}=-\frac{q}{k}---\left(10b\right)$

式中，q为热源强度，k为常数，r表示与热源的距离，θ为角度。

边界条件为：

u(a,θ)＝u₁(1+cos2θ)(11b)

u(b,θ)＝u₂(1+cos2θ)(12b)

其中，a为圆筒内边界，b为圆筒外边界，θ为角度。若采用解析方法则得到的方程的解为：

$\begin{array}{c}u\left(r,\mathrm{\&theta;}\right)=\&lsqb;\frac{u_1\ln a-u_2\ln b}{\ln a-\ln b}+\frac{q_0}{32k}\frac{b_4\ln a-a^4\ln b}{\ln a-\ln b}\&rsqb;+\&lsqb;\frac{u_2-u_1}{\ln a-\ln b}+\frac{q_0}{32k}\frac{a^4-b^4}{\ln a-\ln b}\&rsqb;\ln r\\ +\&lsqb;\frac{a^2{u}_2-b^2{u}_1}{a^4-b^4}+\frac{q_0}{24k}\frac{a^4+a^2{b}^2+b^4}{a^2+b^2}\&rsqb;r^2\cos 2\mathrm{\&theta;}+\&lsqb;\frac{a^2{b}^2\left(a^2{u}_1-b^2{u}_2\right)}{a^4-b^4}-\frac{q_0}{24k}\frac{a^4{b}^4}{a^2+b^2}\&rsqb;\frac{\cos 2\mathrm{\&theta;}}{r^2}\\ +\frac{q_0}{8k}{r}^4\left(\frac{1}{4}+\frac{1}{3}\cos 2\mathrm{\&theta;}\right)\end{array}---\left(13b\right)$

式中，a为圆筒内边界，b为圆筒外边界，θ为角度，q₀、k为常数，r表示与热源的距离。

对于雅克比迭代法、高斯赛德迭代法没有最佳收敛因子问题，在相同收敛精度条件下重复计算，其迭代次数基本不变；对于超松弛迭代、并行超松弛迭代法，在迭代次数为100，网格划分为100×100条件下，松弛因子ω选取1.72时，算法精度分别是1.284×10^-3和9.130×10^-4；ω选取1.80时，算法精度分别是8.732×10^-4和7.557×10^-4，后者的收敛效果比前者更好，表明最佳因子不能依靠经验取得。图3(a)、(b)和(c)表明适应度函数是一个凸函数，具有唯一极值，且极值对应的松弛因子在1.94附近，这个特点与理论分析是一致的。由图3(d)知，随着种群的进化，种群的适应度在不断增加，但到进化后期，其值基本不变，维持在0.99～1范围内，表明进化是一个择优的过程，此外适应度最大且稳定不变也是遗传计算收敛的表现。图3(e)，(f)分别是进化过程中每代的最佳因子与最大迭代误差的变化曲线，可见随着种群的进化，误差逐渐减小，充分体现了“适者生存”的遗传思想，也符合进化的原理。应用本发明经过10代进化，GA寻找最佳因子的过程用时T_GA＝22.971秒，找到的最佳因子ω₀为1.93557977943888，对应的适应度为0.999914589741900。表1、表2表面松弛因子对泊松方程求解的速度与精度影响显著；改进算法能减少迭代次数，节省计算时间，加快方程的求解；算法适合于求解计算量较大、精度要求较高的时域有限差分方程，而且精度要求越高，算法的性能越好，节省的时间也越多。

表1各种算法迭代量

表2各种算法计算时间

Claims (2)

1.一种基于遗传计算的二维泊松方程求解优化方法，其特征在于，包括以下步骤：

（1）采用遗传计算对松弛因子进行全局寻优，适应度函数建模为与迭代次数以及收敛精度有关的多目标适应度函数，式中，表示电场、磁场或温度场点处的位函数，表示当前迭代数；

（2）初始化种群，采用截断选择法与稳态繁殖法相结合对种群进行优胜劣汰的筛选，只保留精英个体，提高种群的多样性；

（3）新个体由父个体的线性插值及非均匀变异产生，交叉概率和变异概率根据自适应遗传算法进行计算；

（4）判断收敛性，遗传算法的收敛条件是迭代次数超过300或者最大适应度连续3代变化都小于，算法收敛时适应度最大值对应的个体即为最佳松弛因子；

（5）若算法收敛，则选用五台处理能力相同的PC机作为硬件平台，一台作为主机，其余作为从机，主机与从机通信，从机之间互不干扰；将遗传计算得到的最佳松弛因子作为并行超松弛迭代算法的松弛因子，进行并行超松弛迭代计算，实现二维泊松方程的快速计算；

若算法不收敛，则返回步骤（2），继续通过截断选择法与稳态繁殖法相结合对种群进行优胜劣汰的筛选；

所述多目标适应度函数为：

；

；

式中c1、c2为正加权系数，满足c1+c2=1，为迭代次数，表示电场、磁场或温度场中点处的位函数，表示松弛因子，表示当前代的最大误差；表示当前代迭代次数的最大值，为正整数。

2.按照权利要求1所述的基于遗传计算的二维泊松方程求解优化方法，其特征在于：所述步骤（2）中，初始化种群时，设置种群大小，染色体长度为15，进化代数，变异算子，截断阈值。