基于粒子群优化算法的散乱数据光顺拟合

1、基于粒子群优化算法的散乱数据光顺拟合    0引言光顺问题是CAGD中的一个经典问题,在工业部门的设计和制造中占有重要地位。光顺是指微调曲线曲面的形状,在其形状变化不超过给定约束的前提下,使曲线曲面变得更加光滑。这里的约束通常是指在曲线曲面上一批数据点的几何约束。在数学放样和工业设计制造中,原始的外形信息一般是离散的型值点,这些数据点往往带有一定的误差。由这些测量点拟合成曲线或曲面时就会产生波动。很多光顺方法所用到的数据点都假设分布在网格线上,即是有序的,而且所采用的参数化方法通常与数据点的分布情况有关。对于无序的数据用以上方法则难以处理好,而且有序化这

2、些点也是困难和耗时的。对于散乱数据点的光顺拟合问题,许多学者做了卓有成效的研究,文献1对光顺准则中特殊参数条件的情况,就光顺问题的SOR法进行了讨论;文献2和3讨论了用有限元方法的求解光顺拟合问题。以上均是在传统优化方法上讨论散乱数据的光顺拟合问题。本文首次采用群体智能粒子群优化算法,增加了计算的自动程度以及客观性。在优化中将模糊集合理论和粒子群优化算法有机地结合起来,采用了优于传统罚函数法的模糊罚函数法处理几何约束。通过对节点序列内在关联性的分析,提出了适合邻域搜索类算法实施的邻域结构,使得曲面光顺的求解在实质上可以应用演化思想优化。试验结果表明,本方法能够有效地实现散乱数据点的光顺拟合。1

3、散乱数据光顺原理及离散化光顺准则涉及到几何外型的美观,其提法有许多种1-5,其中文献2推广了光顺准则,采用更一般的形式,即对于曲面u = u(x,y),(x,y) R2,光顺准则定义为E(u) = 1(ux)2+2(uy)2+1(2ux2)2+2(2uy2)2dxdy(1)式(1)中不仅含有一阶偏导数项,还含有二阶偏导数项,其中、分别相当于拉伸刚度和弯曲刚度,控制着曲面的拉伸和弯曲变形。考虑二维散乱数据的光顺拟合问题,设f(x,y)C1(),其中=a,b×c,d为有界矩形域。由定义域的等距划分:x1=a+ih, i=0,1,I;h =(b-a)I;yi= c+jq,j=0,

4、1,J;q =(d-c)J,得到关于(x,y)的离散点指标集合N= (i,j) |i=0,1,I;j=0,1,J。假设存在N的子集Kij=(lk,mk|k =1,2,N)和有限实数集Z=z1,z2,zN,对任意(lk,mk)Kij N,有zkZ,使f(xlk,ymk) =zk(k=1,2,N),求二元函数u(x,y),使得minu(x,y)E(u)满足ulk,mk=zk,k =1,2,N。采用差分方法将二次泛函E(u)离散化,可得其离散化模型f=minJj=0I-1i=01(ui+1,j-ui,jh)2+Ii=0J-1j=02(ui,j+1-ui,jq)2+Jj=0I-1i=11(ui+1,j

5、-2ui,j+ui-1,jh2)2+Ij=0J-1j=12(ui,j+1-2ui,j+ui,j-1q2)2(2)满足Ck:ulk,mk=zk,k =1,2,N。此即为有约束的优化问题,其中f为目标函数,Ck为约束条件。2优化算法描述2. 1粒子群优化算法及其特点Kennedy等受鸟群觅食行为的启发,于1995年提出粒子群优化算法(PSO)。与进化算法相比, PSO算法保留了基于种群的全局搜索策略,但其所采用的速度和位移相对应的搜索模型操作简单,避免了复杂的进化操作,在解决复杂函数优化等许多实际问题中取得了成功的应用。PSO算法采用以下公式进行单位迭代次数解的更新vi+1=w×vi+c

6、1×rand()×(pi-xi) +c2×rand()×(gi-xi) (3)xi+1=xi+vi+1(4)其中,w为惯性权重;rand()为均匀分布在(0,1)之间的随机数; c1和c2为学习因子。式(3)中的速度vi= (vi1,vi2,vid)决定粒子在解空间内单位迭代次数的位移,它反映了粒子所代表的解xi=(xi1,xi2,xid)单位迭代次数的变化,如式(4)所示。此外,粒子的速度vi被最大速度vmax所限制,即若vi>vmax。则令vi=vmax。每一次迭代,粒子通过2个极值更新其速度和位置。第1个极值是单个粒子从算法迭代初始到当前迭代

7、搜索所生成的最优解,即个体极值pi= (pi1,pi2,pid);第2个则是粒子所在邻域对应的最优解,即邻域极值gi= (gi1,gi2,gid)。粒子在解空间内不断跟踪个体极值和邻域极值进行搜索,直到满足迭代停止条件,即达到规定的迭代次数或满足规定的误差标准。粒子群优化算法的突出特点是:搜索空间是从众多的初始群体开始,搜索过程可以有效地跳出局部极值点,得到问题的全局最优解的概率大大提高。而且,不需要对待寻优参数进行编码,算法中需要选择的参数少,程序实现非常简洁,并在种群数量、寻优速度等方面具有一定的优势,运算量相对较小。2. 2约束条件的处理粒子群优化算法不能直接处理有约束的问题,解决方法之

8、一是利用罚函数法将有约束问题转化为无约束问题。由于粒子群优化算法与传统方法相比具有随机搜索和多个体操作的优点,在粒子群优化算法中对约束条件的处理与传统方法有较大的差别。同时,由于有惩罚项的存在,其违反约束条件的程度会明显地影响搜索进程,考虑到确定性的数学表达式对粒子群优化算法的进化过程影响不大,因此可以将模糊集合理论和粒子群优化算法集成。由于约束条件被转换到模糊域中后,模糊可行集合中既包含可行点又包含不可行点,所以粒子群优化算法能够同时得到可行点和不可行点的信息。在模糊环境下,约束条件由定义域中的模糊集合G定义,用G表示点满足约束条件的程度。由模糊集合理论,当点在可行域中时,其隶属函数G为1。

9、其它情况下,隶属函数的值在0G1的区间范围内。群体中的每一点都和该点在模糊可行域中的隶属度有关,点离可行域越近,在模糊可行域中的隶属度也应该越高,处于可行域中的点应该具有最高的隶属度,为此采用一种离散型隶属函数。对于极小化问题,设Rk是由离散隶属函数所确定的模糊罚函数Rk=int(1-c(x)KDG(x) >0100G(x) =0(5)本研究取KD =10。式(5)中,Rk总是取为整数。这样做的优点是能够很明显地看出该点所处的区域。在设计模糊罚函数中的惩罚项时,注意先将约束条件正规化。对违反约束条件的点,应在它所对应的目标函数后加上惩罚函数项。根据模糊理论中隶属度函数的概念,设计的惩罚函

10、数项为:假设点xk对第i个约束的违反程度为dki。dki=0,gi(xk) <0gi(xk),otherwise设点xk对M个约束的违反程度中最大值为maxD =maxdk1,dk2,dkm,则有惩罚项Rk。Rk反映的不是约束函数的值,而是群体中的点违背约束的程度。它的作用是充分利用那些已经在可行域中的点和离可行域很近,有潜力经过迭代演化以后进入可行域中的点,使排序算法能够顺利进行。Rk求出以后,采用罚函数来处理约束条件,则适应度函数应作如下调整f=f+RkNk=1(ulk,mk-zk)2(6)其中,f为正规化后的目标函数式,Rk是xk点处的惩罚因子值。粒子群优化算法在搜索进化过程中,一

11、般不需要其它外部信息,仅用适应度函数来评价个体或解的优劣,并作为以后进化操作的依据。用粒子群优化算法进行优化的过程也就是对评价函数求极值的过程。2. 3基于邻域结构的搜索简单的随机遍历性搜索,很难达到最优解,采取先对状态序列子列进行局部优化的方法。对状态序列某个子列的局部优化,实质上是在整体关联性不受破坏的前提下对整个状态序列的优化。在这里,局部优化思想引入邻域结构的概念,将局部优化视做某个邻域结构中寻优的过程。更进一步,本文所定义的邻域结构是随机并具有Markov各态遍历性的,因而可发展为真正有效的邻域搜索类算法。在式(6)中,评价函数的各项都是半正定的,泛函指标会随节点间弧长的增长而增加。

12、显然,短的弧长更可能得到较小的泛函指标。当曲面复杂时,很难对各个节点坐标位置有一个整体的了解,那么线性地连接相邻状态的节点正好是最优或者接近最优的可能性较大。用概率分布来描述这一个事实如下。定义若A | sn= sn-1+sn+1-sn-12+e=sn+1+sn-1-sn+12+e, snA (7)则称二元集sn,为状态转移sn-1sn+1的一阶邻域结构。式(7)中e为单位向量,e的方向与sn-1sn+1正交,纯量服从N(0,2)分布,=sn-1sn+1。从定义可见,一阶邻域结构可以遍历一个与sn-1sn+1正交的直线,正态分布P满足在可行域上的Markov各态遍历性。2. 4算法步骤散乱数据

13、点光顺拟合问题的粒子群优化算法流程为:1)初始化所有粒子(群体规模为m)。迭代初值取为ui,j=zk, (i,j) = (lk,mk)Kijrand(),(i,j)N/Kij(8)每个粒子的pi设为初始位置,pi中的最好值设为gi。2)在初值状态序列中以等概率选取节点序列sn,以公式(7)的方式确定其邻域结构sn,并由此产生一个新的sn,形成新的状态序列。3)判断节点序列中每个节点在可行域中的位置,从而得出隶属度确定惩罚因子,将约束问题转化为无约束问题。4)评价每个粒子的适应值。根据目标函数(6)计算每个粒子的适应值.如果优于pi,则gi被当前位置替换。如果所有粒子的gi中有优于gi的,则重新

14、设置gi的索引号。5)根据式(3)、式(4)调整粒子的位置和速度。6)检查终止条件。如果达到最大迭代次数或者最好解停滞不再变化,就终止迭代,否则以当前结果作为初值状态,回到步骤2。3算例将平面矩形区域=77, 194×-81,143剖分为14×15网格,分划为xi=77+ih, i=0,1,13;h=(194-77)/13=9yi=-81+iq, i=0,1,14;q=(143+81)/14=16记ui,j=u(xi,yj), (i,j)N=(i,j) |i=0,表1散乱数据点坐标数据Tab.13-D scattered points before fitting编号k横坐标编号Ik纵坐标编号mk竖坐标值zk1 5 1 3.02 11 2 4.03 7 3 3.54 10 4 3.05 1 5 1.06 8 6 4.07 12 7 6.08 6 8 6.09 2 9 3.010 11 10 9.011 4 11 3.012 13 12 8.013 9 13 8.014 3 14 1.51,2,13;j=0,1,2,14,当(lk,mk)N时,ulk,mk=zk(k =1,2,14)的数据由表1给出。在初值状态中以相等概率