第 23卷第 5期 　　　 2006年 5月 　　 计算机应用与软件 Computer App lications and Software Vol 23, No. 5 May 2006 基于粒子群算法的非线性方程组求解 (电子科技大学软件学院 　四川 成都 610054) 　 (泸州职业技术学院 　四川 泸州 646005) 陈长忆 　叶永春 摘 　要 　　将非线性方程组的求解问题转化为无约束极大极小优化问题 ,并应用一种新的进化计算 ( EC)方法 ———粒子群算法 ( PSO )求解此优化问题 。数值实验的结果验证了该方法的可行性和有效性 。 关键词 　　粒子群算法 　非线性方程组 　极大极小问题 SOL V ING NO NL INEAR SY STEM S O F EQUAT IO NS BASED O N PART ICL E SW ARM O PT IM IZAT IO N Chen Changyi　Ye Yongchun ( S chool of Softw a re, U niversity of E lectron ic S cience & Technology of China, Chengdu S ichuan 610054, Ch ina) (L uzhou V ocationa l & Techn ical College, L uzhou S ichuan 646005, Ch ina) Abstract　　The solving p roblem for nonlinear system s of equations is transformed to the unconstrained m inimax op tim ization p roblem, and then, the particle swarm op tim ization ( PSO) algorithm which is a new evolutionary computation ( EC) method is used for solving the m inimax p roblem. Numerical results show the feasibility and efficiency of this method. Keywords　　Particle swarm op tim ization　Nonlinear system of equation　M inimax p roblem 1　引 　言 众所周知 ,寻求求解非线性方程组的有效方法一直是工程 技术界和应用数学界的难点之一 。近 20年来 ,国内外的许多数 学家对非线性方程组的求解问题进行了大量研究 , 提出了许多 有效解法 , 如消元法 、数值迭代法和进化计算方法等 [1, 2 ]。其 中 ,数值迭代法是最经典和最常用的方法之一 ,但所有的数值迭 Raphson迭代法 (NR 法 )为基础 [ 2 ] 。NR 法具 代法均以 Newton 有局部二阶收敛性 ,且敛速较快 。但 NR法尚有一些弱点 [ 1 ] : ① NR法的每步迭代都要计算 Jacobi矩阵并求解线性方程组 ,工作 量较大 ; ②多数情况下 , NR 法对初值的选择较为敏感 ,要给出 收敛到所需要解的初值较困难 ; ③当 Jacobi矩阵奇异或接近奇 异时 , NR法的计算将无法进行 。 为克服 NR法的以上弱点 ,学者们提出了很多变型算法 。 文献 [ 3 ]提出了求解非线性方程组的极大熵法 ,其核心思想是 借助极大熵原理构造可微熵函数来逼近目标函数 。尽管该方法 的收敛性较好 ,但仍需应用数值迭代法求解 ,其收敛结果依赖于 初值 ,若熵函数的参数选择不当还容易造成上溢 。 近年来 ,进化计算被广泛用于优化问题的求解 ,其中遗传算 法 ( GA)作为一种基于随机性的非数值并行算法得到了更多的 关注 。当前 ,在进化计算领域内又有一种新颖的 、简便的智能算 法 ———粒子群算法 ( Particle Swarm Op tim ization, PSO )受到研究 者的重视 。 PSO是由 J. Kennedy和 R. C. Eberhart[ 4, 5 ]于 1995年 提出的一种新的仿生优化方法 ,其思想最初来源于对鸟群觅食 过程的模拟 ,从这一简单社会系统的行为模型中得到启示 ,由此 建立算法的开发基础 ,并最终发展成为一种有效的优化工具 。 PSO也是基于群体的非数值并行算法 ,整个过程毋须导数信息 。 同传统的优化方法相比 , PSO具有更强的全局优化能力 ,也能较 快地收敛于可接受解 [ 6, 7 ] 。凡能用 GA 求解的问题均能用 PSO 求解 ,而且 PSO不需编码和译码 ,亦不需要诸如选择 、交叉 、变 异和保留等遗传操作 ,因此控制参数少 ,算法简单 ,效率较高 。 2　粒子群算法 ( PSO ) 假设在一个群体中有 pop size个粒子 ,每个粒子是 n维空间 R n 中的一个个体 , 不同的个体具有不同的位置 X j ( j = 1, 2, …, pop size) ,不同的位置 X j 对应于不同的与优化目标函数值相关 的个体适应度函数值 fitness ( X j ) 。 对于第 j个粒子 , 以 X j [ xj1 , xj2 , …, xjn ]T 表示其位置 , 它在 空间中以某个速度 V j [ vj1 , vj2 , …, vjn ]T 飞行 ,它所经历的最好位 置为 Pj = [ pj1 , pj2 , …, pjn ]T ,也称为 Pbest。群体中所有粒子所经 历的最好位置为 Pg = [ pg1 , pg2 , …, pgn ]T ,也称为 gbest。若已知第 k代第 j个粒子的速度 V ( k) ,则第 ( k + 1)代第 j个粒 子的速度及位置为 : V ( k + 1) …, pop size) ) 　 ( j = 1, 2, ( 1) ( 2) 式中 : w ( k) ———惯性权重系数 , 为迭代次数的函数 , 且随迭代次 数线性减少 ,即 [ 6 ] : ) + c2 rj2 ( Pg - X ( k) j 及位置 X ( k) j + c1 rj1 ( Pj - X ( k) j =w ( k) V ( k) X ( k + 1) j j = X ( k) j +V ( k + 1) j j j ( 3) 这里 wmax, wm in ———初始 、终止惯性权重 ; kmax ———最大迭代代数 ; w ( k) =wmax - k (wmax - wm in ) / kmax 收稿日期 : 2005 - 05 - 18。陈长忆 ,讲师 ,主研领域 :计算智能 ,算法 设计与分析 。 

　 138　　 　　　 计算机应用与软件 2006年 c1 ———粒子自身加速度权重系数 ,一般在 0～2之间取值 ; c2 ———全局加速度权重系数 ,一般在 0～2之间取值 ; rj1 , rj2 ———[ 0, 1 ]范围内两个相互独立的 、均匀分布的随机数 ,即 rj1 , rj2 ～U ( 0, 1) 。 为了减少在搜索过程中粒子离开搜索空间的可能性 , V 通 常被限定于一定的范围内 ,即 - Vmax ≤V ≤Vmax。当然 ,亦可根据 具体情况 ,将搜索空间限定于一定的范围内 ,即 Xm in ≤X ≤Xmax。 PSO的计算流程如下 : Step1 初始化一群粒子 (种群规模为 pop size) , 包括随机位 置和速度 , k = 0; Step2 计算每个粒子的适应度值 fitness ( X ( k) ) ; Step3 对每个粒子 ,将其适应度值 fitness ( X ( k) )与其历史最 好位置 Pj的适应度值 fitness ( Pj )比较 ,若较好 ,则将其作为当前 的最好位置 Pj; j j Step4 对每个粒子 ,将其适应度值 fitness ( X ( k) )与全局所经 历的最好位置 Pg 的适应度值 fitness ( Pg )比较 , 若较好 , 则重新 设置 Pg; j Step5 用式 ( 1)和式 ( 2)计算 V ( k + 1) j 和 X ( k + 1) j 其中 : ( j = 1, 2, …, popsize) ; Step6 若达到终止条件 (通常为足够好的适应值或达到预 设最大代数 kmax ) , 则结束 , 返回当前全局最优个体 Pg 即为结 果 ;否则 , k = k + 1,转 Step2。 3　求解非线性方程组的极大极小模型 非线性方程组的一般数学模型为 : ( 4) 式中 : X = [ x1 , x2 , …, xn ]T ∈R n为待求的 n个未知量 ; f ( X ) = [ f1 ( X ) , f2 ( X ) , …, fm ( X ) ]T 为变量 X 的 m 维向量值函数 。 f ( X ) = 0 应用无约束优化方法求解问题 ( 4)时 , 通常将其化为非线 性最小二乘问题 : m X∈R n ∑ m in fi 若以 ‖·‖p 表示向量空间 p i =1 2 ( X ) ( 5) 范数 ,则问题 ( 4)转换为求解 : f ( X ) p m in X∈R n 当 p = 2时 ,即为问题 ( 5) ;当 p = ∞时 ,问题 ( 4)转化为无约 束极大极小优化模型 : m in X∈R n m ax 1≤i≤m { | fi ( X ) | } ( 6) 问题 ( 6)为非光滑优化问题 ,应用 PSO求解非常方便 ,不必 构造可微光滑函数 (如文献 [ 3 ]中的极大熵函数 )来逼近目标函 数 。因此 ,应用 PSO求解非线性方程组式 ( 4)时 ,以式 ( 6)作为 粒子的适应度评价函数 ,且函数值愈小表示粒子的适应度愈高 。 基于优化模型 (6) ,笔者编写了应用 PSO 求解非线性方程组的 C程序 NLSEPSO. C。 4　数值实验 例 1[3 ] 　 1 - x2 + 1 = 0 f1 ( X ) = x2 f2 ( X ) = x1 - cos[ (π /2) x2 ] = 0 应用 PSO求解该问题 ,调用程序 NLSEPSO. C, 取向量维数 n = 2,群体规模 popsize = 30,收敛条件 ε = 1 ×10 - 6 ,最大迭代代 数 kmax = 500, 初始惯性权重 wmax = 1. 2, 终止惯性权重 wm in = 0. 1,加速度权重 c1 = c2 = 1. 8,最大速度 vmax = 0. 1,独立运行 20 次 (全部收敛 )得到全部 3组不同实数解 : X (1) = [ 0. 000 000, 1. 000 000 ]T , X (2) = [ - 0. 707 107, 1. 500 000 ]T , X (3) = [ - 1. 000 000, 2. 000 000 ]T (文献 [ 3 ]未得到第 3组解 ) 。 文献 [ 3 ]指出 ,用 NR 法的连续形式求解该问题时 ,其收敛 性不稳定 ;而用 Euler有限差分法时 ,该问题为 Stiff问题 。应用 本文给出的方法很容易地得到了全部实数解且收敛性稳定 。 例 2[2 ] 　 f1 ( X ) = x1 f2 ( X ) = x2 f2 ( X ) = x1 2 x2 2 x3 2 x3 2 + x1 3 + x2 2 + x3 2 + x2 - 10 = 0 2 + x3 - 20 = 0 2 + x1 - 30 = 0 本文在复数域内求解此方程组 ,将所有复数设为 xd = xR d + jxId的形式 ( xR d , xId为 xd的实部 、虚部 。 d = 1, 2, 3, 为向量分量的 下标 ) ,再将之代入原方程组 ,并令实部 、虚部分别相等 ,即可得 到与原方程组等价的新方程组 (未知量 、方程数目均为原来的 2 倍 ) 。应用 PSO 求解该问题 , 调用程序 NLSEPSO. C, 取 n = 6, popsize = 200,ε= 1 ×10 - 5 , kmax = 1 000, wmax = 1. 2, wm in = 0. 1, c1 = c2 = 1. 8, vmax = 0. 1,独立计算 100次 (全部收敛 )得到全部 16 组不同解 。数值结果如表 1所示 。 表 1　例 2的数值结果 No. xR1 x11 xR2 x12 xR3 x13 1 2 3 4 5 6 7 8 9 10 11 12 13 14 1. 745 226 0. 000 000 - 1. 684 075 0. 000 000 - 2. 642 670 0. 000 000 1. 676 546 0. 000 000 1. 431 257 0. 000 000 2. 726 243 0. 000 000 1. 559 926 0. 000 000 1. 569 841 0. 000 000 - 2. 878 096 0. 000 000 1. 849 814 0. 000 000 - 1. 540 316 0. 000 000 2. 523 135 0. 000 000 - 1. 785 181 0. 000 000 - 1. 627 434 0. 000 000 - 2. 755 291 0. 000 000 - 1. 897 347 0. 000 000 - 1. 479 462 0. 000 000 2. 633 310 0. 000 000 - 1. 721 720 0. 000 000 1. 381 135 0. 000 000 2. 828 746 0. 000 000 - 1. 597 308 0. 000 000 1. 523 864 0. 000 000 - 2. 982 811 0. 000 000 0. 055 282 2. 214 733 0. 130 669 - 1. 734 292 - 0. 015 455 - 2. 771 252 0. 055 282 - 2. 214 733 0. 130 669 1. 734 292 - 0. 015 455 2. 771 252 0. 178 236 2. 216 527 0. 007 389 1. 746 822 0. 173 873 - 2. 742 207 0. 178 236 - 2. 216 527 0. 007 389 - 1. 746 822 0. 173 873 2. 742 207 0. 033 662 2. 220 887 0. 083 794 1. 739 762 0. 049 691 2. 763 258 0. 033 662 - 2. 220887 0. 083 794 - 1. 739 762 0. 049 691 - 2. 763 258 15 - 0. 089 188 2. 229 362 0. 054 218 - 1. 738 054 0. 239 162 2. 742 226 16 - 0. 089 188 - 2. 229 362 0. 054 218 1. 738 054 0. 239 162 - 2. 742 226 5　结 　论 将非线性方程组的求解问题转化为无约束极大极小优化问 题 ,并应用一种新的进化计算方法 ———粒子群算法求解此优化 问题 。由于该方法是基于随机性和群体的非数值并行算法 ,其 收敛结果不依赖于初值 ,且整个过程毋须导数信息 ,因此收敛性 较好 。对于多解的非线性方程组 ,可以通过独立计算多次来得 到多个 /完全数值解 。数值实验表明 , 应用 PSO 求解非线性方 程组具有算法简单 、控制参数少且能求出多个 /完全解等优点 , 是一种值得推广的求解非线性系统的实用算法 。 

2 2 2 2 2 　 第 5期 　　　 陈长忆等 :基于粒子群算法的非线性方程组求解 139　　 参 考 文 献 [ 1 ] 蔡大用、白峰杉 ,现代科学计算 [M ] ,北京 :科学出版社 , 2001. [ 2 ] 张纪元 ,机械学的数学方法 [M ] ,上海 :上海交通大学出版社 , 2003. [ 3 ] 孔敏 、沈祖和 ,“解非线性方程组的极大熵方法 [ J ] ”,《高等学校计 算数学学报 》, 1999, 21 (1) : 1～7. [ 4 ] Kennedy J, Eberhart R. Particle swarm op tim ization [ A ]. Proceedings of IEEE International Conference on Neural Networks ( ICNN ’95 ) [ C ]. Perth,WA, Australia, 1995, 1942～1948. [ 5 ] Eberhart R, Kennedy J. A new op tim izer using particle swarm theory [A ]. Proceedings of Sixth International Symposium on M icro Machine and Human Science [ C ]. Nagoya, Japan: Nagoya Municipal Industrial Research Institute, 1995, 39～43. [ 6 ] Shi Yu hui, Eberhart R. Parameter selection in particle swarm op tim iza tion[ A ]. Proceedings of 7 th Annual Conference on Evolutionary Pro gramm ing[ C ]. March 1998, 591～601. [ 7 ] Ray T, L iew K M. A swarm with an effective information sharing mecha nism for unconstrained and constrained single objective op tim ization p roblem s[A ]. Proceedings of IEEE International Conference on Evolu tionary Computation[ C ]. South Korea: IEEE Press, Seoul, 2001, 75 ～ 80. (上接第 27页 ) 4 3　分析与死锁避免 将上面对哲学家就餐问题的 FSP描述输入标号迁移系统 分析器 ( labelled transition system analyzer, LTSA ) ,先进行安全性 检查 ,分析器的输出如图 2所示 。 图 2　安全性检测输出 接着进行活动性检查 ,分析器的输出如图 3所示 。 图 3　活动性检测输出 可以看到 , LTSA很快就检测出这种算法中存在死锁和无法 继续进行的活动 ,因此该算法不满足安全性和活动性要求 ,同时 LTSA也给出了发生死锁的一种情形 ,即 5个哲学家同时拿起左 边的一支叉子 ( left. get) ,然后等待右边的叉子 ,因此谁也等不 到他们需要的叉子 ,即发生了死锁 。 因此 ,为避免发生死锁 ,我们必须修改 FSP描述 ,对哲学家 进行编号 ,使编号为偶数的哲学家先拿起他们左边的叉子 ,而编 号为奇数的哲学家则先拿起他们右边的叉子 ,以消除循环等待 的可能 (当然还有其它的方法来避免死锁 ) 。修改后的 FSP描 述如表 4所示 。 表 4　修改后的哲学家就餐问题描述 PH IL ( I = 0) = (when ( I% 2 = = 0) sitdown - > left. get - > right. get - > eat - > left. put - > right. put - > arise - > PH IL | when ( I% 2 = = 1 ) sitdown - > right. get - > left. get - > eat - > left. put - > right. put - > arise - > PH IL) FORK = ( get - > put - > FORK) | |D INERS (N = 5) = forall[ i: 0. . N - 1 ] (phil[ i]: PH IL ( i) | | {phil[ i]. left, phil[ ( ( i - 1) +N) %N ]. right}: : FORK) 现在对修改后的描述进行安全性检查 ,分析器的部分输出 如下 : 　　No deadlocks/ errors 　　Analysed in 157m s 进行活动性检查 ,分析器的部分输出如下 : 　　No p rogress violations detected. 　　Process checked in: 266m s. 可见 ,修改后的算法已经满足安全性和活动性的要求 。 5　结束语 在并发和分布式系统的设计过程中 ,必须满足安全性和活 动性这两个基本性质 。但是由于状态空间和可能发生的迁移数 非常巨大 ,因此手工的分析是困难的 。本文介绍的基本 FSP模 型的描述方法及其分析工具 LTSA能非常迅速地检测并发系统 和分布式系统的安全性和活动性 ,使人们在设计这些系统或算 法的过程中能够尽早发现和修正结构上的问题和错误 ,具有重 要的实际意义 。分析工具 LTSA 采用 JAVA 编写 ,可以从 http: ∥wwwdse . doc. ic. ac. uk /～ jnm 获得 ,也可以自行开发以实现 更强的功能 。 参 考 文 献 [ 1 ] Schneider F. B. , Andrews G. R. , Concep ts for Concurrent Programm ing [ J ] , LNCS224, Sp ringer2Verlag, 1985. 669～716. [ 2 ] Gerard Tel. , Introduction to D istributed A lgorithm s ( second edition) [M ] , Cambridge University Press, Cambridge, U. K. , 2000. [ 3 ] Cheung S. C. and Kramer J. , Checking Subsystem Safety Properties in Compositional Reachability Analysis, 18 th IEEE Int. Conf. on Software Engineering( ICSE 18) , Berlin, 1996, 144～154. [ 4 ] Cheung S. C. , Giannakopoulou D. , and Kramer J. , Verification of L ive ness Properties using Compositional Reachability Analysis, 6 th European Software Engineering Conference /5 th ACM SIGSOFT Symposium on the Foundations of Software Engineering ( ESEC / FSE 97 ) , Zurich, Sep t. 1997, LNCS 1301, ( Sp ringer Verlag) , 1997, 227～243. [ 5 ] Jeff Magee and Jeff Kramer, Concurrency State Models & Java Program s [M ]. W iley, 1999. [ 6 ] ANDREW S TANENBAUM, ALBERTSWoodhull, Operating system s: de sign and imp lementation ( second edition) [M ] , Prentice Hall, 1997.