查看: 28012|回复: 129|关注: 33

[我分享] 基于全局优化工具箱求解非线性方程组数值解(无需初值)

  [复制链接]

论坛优秀回答者

权威

3630 麦片

财富积分



20

主题

3779

帖子

769

最佳答案
  • 关注者: 431
发表于 2017-2-26 19:59:55 | 显示全部楼层 |阅读模式
本帖最后由 wuyou136 于 2017-2-26 20:24 编辑

通常我们求解非线性方程组都是用 fsolve(Optimization Toolbox),但是fsolve 有两大局限性:
1、需要自己判断搜索初值。fsolve 属于局部优化函数,对于初值非常敏感,在对目标函数性质不明确的时候使用起来非常不方便,可能自己尝试很多次初值都不能得到收敛结果。
2、不能约束解的范围。当目标函数有多解时,如果我们需要对 fsolve 限定求解范围,通常是将搜索初值取在求解范围之内,不过这有时候也可能导致最终求得的解并不在我们想要的范围之内(初值与解的范围没有必然关系)。

这里利用 Global Optimization Toolbox 中的 particleswarm+GlobalSearch,实现非线性方程组的全局搜索,可以解决上述两个问题。
其原理是先利用 particleswarm 得到一个接近于解的搜索初值,然后将目标函数限定为 fmincon 的非线性约束,最后利用 GlobaSearch 寻找目标解。
以 :


这个方程组为例,它一共有4组解:
x = (–1,–2)
x = (10,–2),
x = (–1,20),
x = (10,20).
如果我们想要得到 x=(-1,-2) 这组解:
  1. function main
  2. lb = [-inf -inf]; ub = [0 0]; %限定求解范围
  3. x0 = particleswarm(@(x)sum(myfun(x).^2),2,lb,ub); %利用 particleswarm 得到一个搜索初值
  4. opts = optimoptions(@fmincon,'Algorithm','interior-point','Display','off');
  5. problem  = createOptimProblem('fmincon','x0',x0,'objective',@(x)0,'lb',lb,'ub',ub,'nonlcon',@(x)constrain(x,@myfun));
  6. gs = GlobalSearch;
  7. [x,~,exitflag] = gs.run(problem) %全局搜索


  8. function F = myfun(x)

  9. F(1) = (x(1)+1)*(10-x(1))*(1+x(2)^2)/(1+x(2)^2+x(2));
  10. F(2) = (x(2)+2)*(20-x(2))*(1+x(1)^2)/(1+x(1)^2+x(1));

  11. % fmincon 的非线性约束
  12. function [c ceq] = constrain(x,f)
  13. ceq = f(x);
  14. c =[];
复制代码
  1. x =
  2. -1.0000   -2.0000
复制代码
附件将上述过程封装成一个函数,使用方法:

[x,fval,exitFlag] = GlobalSolve(@myfun,n)

[x,fval,exitFlag] = GlobalSolve(@myfun,n,lb,ub)
myfun 为目标函数,n 为求解变量的个数,lb 和 ub 为解的限定范围,默认为 -inf 到 inf

目前这种方法的局限性:
1、不能求解复数方程。
2、如果目标函数不连续,可能会导致求解失败。

大家如果有好的意见,欢迎提出讨论,谢谢!

GlobalSolve.m

898 Bytes, 下载次数: 879

论坛优秀回答者

权威

3630 麦片

财富积分



20

主题

3779

帖子

769

最佳答案
  • 关注者: 431
 楼主| 发表于 2017-2-26 20:11:19 | 显示全部楼层
本帖最后由 wuyou136 于 2017-2-26 20:14 编辑

需要直接使用的看二楼:
将 GlobalSolve 下载放到搜索路径之后,输入:
  1. F = @(x)[(x(1)+1)*(10-x(1))*(1+x(2)^2)/(1+x(2)^2+x(2));
  2. (x(2)+2)*(20-x(2))*(1+x(1)^2)/(1+x(1)^2+x(1));] %目标函数
  3. [x,fval,exitFlag] = GlobalSolve(F,2)
复制代码
修改你的目标函数即可


论坛优秀回答者

权威

3630 麦片

财富积分



20

主题

3779

帖子

769

最佳答案
  • 关注者: 431
 楼主| 发表于 2017-2-27 17:53:52 | 显示全部楼层
本帖最后由 wuyou136 于 2017-2-28 13:00 编辑

更新:增加多解求解功能。以1L的方程组为例:

  1. F = @(x)[(x(1)+1)*(10-x(1))*(1+x(2)^2)/(1+x(2)^2+x(2));
  2. (x(2)+2)*(20-x(2))*(1+x(1)^2)/(1+x(1)^2+x(1));] %目标函数
  3. [x,fval,exitFlag,multiSol] = GlobalSolve(F,2)
复制代码


第4个返回值即求得的多解:
multiSol{1} =

   -1.0000   -2.0000



multiSol{2} =

   10.0000   20.0000



multiSol{3} =

    -1    20



multiSol{4} =

   10.0000   -2.0000




GlobalSolve.m

1.85 KB, 下载次数: 322

论坛优秀回答者

权威

3630 麦片

财富积分



20

主题

3779

帖子

769

最佳答案
  • 关注者: 431
 楼主| 发表于 2017-3-11 17:49:25 | 显示全部楼层
本帖最后由 wuyou136 于 2017-3-11 17:50 编辑

更新:
1、增加第5个输入参数。第5个输入参数含义为 GlobalSolve 的求解能力,默认值是100。这个值给的越高,求解能力越强,但是耗费的时间也越长。

2、可以求复数解。

3、第5个参数给定相同的情况下,求解时间缩短,效率提升大约一倍。

GlobalSolve1.m

1.85 KB, 下载次数: 504

新手

10 麦片

财富积分


050


0

主题

2

帖子

0

最佳答案
发表于 2017-3-11 10:52:48 | 显示全部楼层
楼主真是好人,这种方法很不错。

论坛优秀回答者

退役版主

2206 麦片

财富积分



46

主题

2115

帖子

303

最佳答案
  • 关注者: 160
发表于 2017-3-11 17:16:18 | 显示全部楼层
可以试一下特殊数学函数的多个根,比如贝塞尔函数等。
请点击"回复此楼",否则我将无法收到回帖提醒。
问题如果比较复杂或较难,请邮箱联系kimist@qq.com

论坛优秀回答者

权威

3630 麦片

财富积分



20

主题

3779

帖子

769

最佳答案
  • 关注者: 431
 楼主| 发表于 2017-3-11 17:46:02 | 显示全部楼层
kastin 发表于 2017-3-11 17:16
可以试一下特殊数学函数的多个根,比如贝塞尔函数等。

尝试了一下,GlobalSolve 对于一元函数求解效果比较好,基本能找出所有解。比如对于 [0 50] 范围内的贝塞尔函数,一共17个根,GlobalSolve 求解没有漏解:

[~,~,~,ms,mf] = GlobalSolve1(@(x)besselj(0.1,x),1,0,50)
X.jpg

论坛优秀回答者

退役版主

2206 麦片

财富积分



46

主题

2115

帖子

303

最佳答案
  • 关注者: 160
发表于 2017-3-11 18:45:03 | 显示全部楼层
wuyou136 发表于 2017-3-11 17:46
尝试了一下,GlobalSolve 对于一元函数求解效果比较好,基本能找出所有解。比如对于 [0 50] 范围内的贝塞 ...

漏掉倒不大可能,我想的是不知道精度如何(Mathematica有专门的零点函数的作为标准对比)。
另外,可以选择一些高振荡的函数来测试一下鲁棒性。
请点击"回复此楼",否则我将无法收到回帖提醒。
问题如果比较复杂或较难,请邮箱联系kimist@qq.com

论坛优秀回答者

权威

3630 麦片

财富积分



20

主题

3779

帖子

769

最佳答案
  • 关注者: 431
 楼主| 发表于 2017-3-11 21:20:30 | 显示全部楼层
kastin 发表于 2017-3-11 18:45
漏掉倒不大可能,我想的是不知道精度如何(Mathematica有专门的零点函数的作为标准对比)。
另外,可以选 ...

版主,我尝试了下,GlobalSolve 和 fsolve 的求解精度是完全一样的,它的求解精度取决于 lsqnonlin 的求解精度,可以通过参数选项来控制。
以前面贝塞尔函数为例,如果误差参数我都给定为默认值的话:
  1. n = 1; lb = 0; ub = 20; obj_fun = @(x)besselj(0.1,x); start_num = 100;

  2. options = optimoptions('particleswarm','display','off');
  3. x0 = particleswarm(@(x)sum(obj_fun(x).^2),n,lb,ub,options);
  4. opts = optimoptions(@lsqnonlin,'FunctionTolerance',1e-6,'StepTolerance',1e-6,'OptimalityTolerance',1e-6);
  5. problem  = createOptimProblem('lsqnonlin','x0',x0,'objective',obj_fun,'lb',lb,'ub',ub,'options',opts);
  6. gs = MultiStart('FunctionTolerance',1e-2,'XTolerance',1e-2);
  7. [x,~,exitflag,~,allmins] = gs.run(problem,start_num);
  8. fval = obj_fun(x);
  9. for ii = 1:numel(allmins)
  10.     multiSol(ii) = allmins(ii).X;
  11. end
  12. multiSol = sort(multiSol);
  13. for ii = 1:numel(allmins)
  14.     multiFval(ii) = obj_fun(multiSol(ii));
  15. end
  16. multiSol
  17. multiFval
复制代码

multiFval =

                   0   0.000000000001342  -0.000000000000106   0.000000000057395  -0.000000000000467  -0.000000000001010  -0.000000000041795   0.174759635771728
修改一下误差参数:
  1. n = 1; lb = 0; ub = 20; obj_fun = @(x)besselj(0.1,x); start_num = 100;

  2. options = optimoptions('particleswarm','display','off');
  3. x0 = particleswarm(@(x)sum(obj_fun(x).^2),n,lb,ub,options);
  4. opts = optimoptions(@lsqnonlin,'FunctionTolerance',1e-30,'StepTolerance',1e-30,'OptimalityTolerance',1e-30);
  5. problem  = createOptimProblem('lsqnonlin','x0',x0,'objective',obj_fun,'lb',lb,'ub',ub,'options',opts);
  6. gs = MultiStart('FunctionTolerance',1e-2,'XTolerance',1e-2);
  7. [x,~,exitflag,~,allmins] = gs.run(problem,start_num);
  8. fval = obj_fun(x);
  9. for ii = 1:numel(allmins)
  10.     multiSol(ii) = allmins(ii).X;
  11. end
  12. multiSol = sort(multiSol);
  13. for ii = 1:numel(allmins)
  14.     multiFval(ii) = obj_fun(multiSol(ii));
  15. end
  16. multiSol
  17. multiFval
复制代码

multiFval =

   1.0e-15 *

                   0   0.089815815890432  -0.029356060856082  -0.090921069177840  -0.073607068724981   0.298380430774017   0.162834000060406
另外,我把 GlobalSolve 得到的解作为 fsolve 的初值,选项参数都给一样的话,两者得到的解是完全一致的,所以说两者的求解精度完全一致。不知道 Mathematica 是如何以专门的零点函数作为对比的?我尝试找了下,没找到你说的零点函数

对于高震荡函数的话,不知道你说的是不是附图这种很多零点情况?还是说高震荡函数但是只有个别零点,GlobalSolve 是否能收敛到这些零点?
B.jpg

论坛优秀回答者

退役版主

2206 麦片

财富积分



46

主题

2115

帖子

303

最佳答案
  • 关注者: 160
发表于 2017-3-13 15:22:43 | 显示全部楼层
wuyou136 发表于 2017-3-11 21:20
版主,我尝试了下,GlobalSolve 和 fsolve 的求解精度是完全一样的,它的求解精度取决于 lsqnonlin 的求 ...

这几天比较忙,没来得及回复,见谅。

一般我用Mathematica用于精度验证,因为它相当于matlab的符号函数版,但是速度更快。比如matlab计算调和级数,如果不用符号函数,直接数值计算求和相加到很多项之后,累计误差会很大,但Mathematica不会。对于一些数学上的特殊函数,更是算法得到最大的优化,精度上有保证。

Mathematica中有如下函数
BesselJZero[n, k] 表示贝塞尔函数Jn(x)的第k个根。
BesselJZero[n, k,x0] 表示贝塞尔函数Jn(x)大于x0的第k个根。
当然,k如果是数组,它会自动给出向量化的答案。

对于高振荡函数,可以是逐渐振荡加大的特性,这样就能很好对于算法收敛速度以及收敛精度进行评估,从而选择更好的算法或者选项来求解。选择一些有解析根的例子来检验,比如cos(20 x^2)这种。

因为我的版本比较低,所以没法这里给出检验了……
请点击"回复此楼",否则我将无法收到回帖提醒。
问题如果比较复杂或较难,请邮箱联系kimist@qq.com
您需要登录后才可以回帖 登录 | 注册

本版积分规则

关闭

站长推荐上一条 /3 下一条

快速回复 返回顶部 返回列表