查看: 36650|回复: 148|关注: 41

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

  [复制链接]

论坛优秀回答者

权威

3660 麦片

财富积分



20

主题

3779

帖子

770

最佳答案
  • 关注者: 465
发表于 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, 下载次数: 1077

回复主题 已获打赏: 0 积分

举报

论坛优秀回答者

权威

3660 麦片

财富积分



20

主题

3779

帖子

770

最佳答案
  • 关注者: 465
 楼主| 发表于 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)
复制代码
修改你的目标函数即可


回复此楼 已获打赏: 0 积分

举报

论坛优秀回答者

权威

3660 麦片

财富积分



20

主题

3779

帖子

770

最佳答案
  • 关注者: 465
 楼主| 发表于 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, 下载次数: 425

回复此楼 已获打赏: 0 积分

举报

论坛优秀回答者

权威

3660 麦片

财富积分



20

主题

3779

帖子

770

最佳答案
  • 关注者: 465
 楼主| 发表于 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, 下载次数: 642

回复此楼 已获打赏: 0 积分

举报

论坛优秀回答者

权威

3660 麦片

财富积分



20

主题

3779

帖子

770

最佳答案
  • 关注者: 465
 楼主| 发表于 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
回复此楼 已获打赏: 0 积分

举报

论坛优秀回答者

权威

3660 麦片

财富积分



20

主题

3779

帖子

770

最佳答案
  • 关注者: 465
 楼主| 发表于 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
回复此楼 已获打赏: 0 积分

举报

论坛优秀回答者

权威

3660 麦片

财富积分



20

主题

3779

帖子

770

最佳答案
  • 关注者: 465
 楼主| 发表于 2017-3-13 20:34:35 | 显示全部楼层
kastin 发表于 2017-3-13 15:22
这几天比较忙,没来得及回复,见谅。

一般我用Mathematica用于精度验证,因为它相当于matlab的符号函数版 ...

谢谢,你在论坛都回答都非常有质量,我也是很希望能得到你的建议。:)  

对于方程求解的精度问题,是否只需将求得的解代回到原方程就可以了,也就是函数返回的 fval,而不需要再与 mathematica 的解析解进行对比?因为 double 型最多精确到小数点后16位左右,所以这种数值方法最高是可以精确到小数点后16位的,比如对于这里的贝塞尔函数,fval 在10^-16量级,已经达到了 double 型数值计算的最高精度,也就是说 GlobalSolve 通过调节误差控制参数是可以达到16位精度的,不知道这样理解是否正确?


另外你举例的高震荡函数来对收敛速度和收敛精度进行评估我还是没有理解,这怎么对方程求解进行评估?是指例如与前100个解析根进行对比吗?

回复此楼 已获打赏: 0 积分

举报

论坛优秀回答者

权威

3660 麦片

财富积分



20

主题

3779

帖子

770

最佳答案
  • 关注者: 465
 楼主| 发表于 2017-3-14 17:09:38 | 显示全部楼层
kastin 发表于 2017-3-13 22:17
谢谢你这么看重~

之前我也曾想过如何求解方程在某个范围内的所有根,但是由于没有一般的全局搜索算法,所 ...

谢谢建议。我尝试计算了一下你给的这个函数,可以通过调节 MultiStart 的初值点数,来增强搜索能力找寻所有的根。

在 n 的值比较小的时候,找寻所有根比较容易。但是一旦 n 的值增加,要想不遗漏根的话,初值点数要给的很大才行,并且它漏掉的根比较有规律,都是最靠近原点的那个根,我觉得这也是 GlobalSolve 后续可以改进的地方。

下面附上我测试的一些数据,对于不同的 n 值,如果要求得所有的根,GlobalSolve 所耗费的时间:

  1. a = 1; n = 10; lb = -n-5; ub = n+5;
  2. syms x
  3. fun = sum(1/x+2*x./(x^2-(1:n).^2))-a;
  4. fun = matlabFunction(fun);
  5. fplot(fun,[lb,ub]); hold on
  6. %%
  7. x = linspace(lb,ub,2e3); y = fun(x);
  8. [~,locs] = findpeaks(-abs(y));
  9. h = plot(x(locs),-0.5,'bo'); %蓝色标记代表 findpeaks 找到的根
  10. %%
  11. tic
  12. [~,~,~,ms,mf] = GlobalSolve(fun,1,lb,ub,100);
  13. toc
  14. plot(ms,zeros(size(ms)),'ro') %红色标记代表 GlobalSolve 找到的根
复制代码


n = 10
时间已过 2.036699 秒。

n = 30
时间已过 5.845213 秒。

n = 70
时间已过 25.255925 秒。

上面都是以默认精度进行计算的,fval 在 10^-5 量级,我把计算精度调至最高,发现 fval 是在 10^-11 量级,相应的计算时间也会增加。

n = 20

n = 20

n = 50

n = 50

n = 100

n = 100
回复此楼 已获打赏: 0 积分

举报

论坛优秀回答者

权威

3660 麦片

财富积分



20

主题

3779

帖子

770

最佳答案
  • 关注者: 465
 楼主| 发表于 2017-3-14 19:10:21 | 显示全部楼层
kastin 发表于 2017-3-14 17:25
目前看来效果还是很不错的,不过时间消耗还是比较大的。

我有几个想法:

缩小 lb 和 ub 可以降低搜索时间,就相当于在单位区间内增加初值点数。我觉得你分析的很对,漏掉靠近原点的根应该是因为 1/x 这一项的原因。我用来计算零点的话,就不会出现漏掉靠近原点的根的情况。之前也是想到缩小搜索区间来重新计算靠近原点的根,刚尝试了下,缩小搜索区间的话,也没有漏掉靠近原点的根了。

要提升 GlobalSolve 搜索效率,很重要的一点就是估算 lb 和 ub,这个要根据实际情况来了,lb 和 ub 给的越精确,GlobalSolve 的搜索效率也就越高。

回复此楼 已获打赏: 0 积分

举报

论坛优秀回答者

权威

3660 麦片

财富积分



20

主题

3779

帖子

770

最佳答案
  • 关注者: 465
 楼主| 发表于 2017-3-16 17:14:04 | 显示全部楼层
本帖最后由 wuyou136 于 2017-3-16 17:21 编辑
kastin 发表于 2017-3-14 21:59
既然这样,如果可以事先确定根的区间(包括人为给定区间),那么能有效提高效率和降低遗漏;对于复杂函数 ...

抱歉回复晚了。还是以前面的震荡曲线为例,现在想将区间从 [-30 30] 由算法缩小到 [-10 10],但是二分法我没想明白怎么做可以将这个区间判断出来。
如果要有什么办法,能够判断出来[-30 -10]这一段函数值与 -30 这个端点函数值符号都一致的话那就好了,可以缩小区间同时也可以应用到多元函数里面。

M.jpg



回复此楼 已获打赏: 0 积分

举报

您需要登录后才可以回帖 登录 | 注册

本版积分规则

关闭

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

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