楼主: wuyou136

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

  [复制链接]

论坛优秀回答者

权威

3630 麦片

财富积分



20

主题

3779

帖子

769

最佳答案
  • 关注者: 430
 楼主| 发表于 2017-3-17 09:18:07 | 显示全部楼层
寒山寺 发表于 2017-3-17 09:16
楼主这个函数怎么使啊,为什么我在用这个函数验证时,出现下面这样的错误呢? ...

这是版本问题,你用的时候把这一句:
  1. gs = MultiStart('FunctionTolerance',1e-2,'XTolerance',1e-2);
复制代码

修改成:
  1. gs = MultiStart();
复制代码

论坛优秀回答者

退役版主

2206 麦片

财富积分



46

主题

2115

帖子

303

最佳答案
  • 关注者: 160
发表于 2017-3-20 11:07:58 | 显示全部楼层
wuyou136 发表于 2017-3-11 17:49
更新:
1、增加第5个输入参数。第5个输入参数含义为 GlobalSolve 的求解能力,默认值是100。这个值给的越高 ...

x,feval指的是第一个解吗?既然如此,为何不直接给multiSol和multiFeval呢?

刚刚试了复数解,不知道为何GlobalSolve(@(x)x^2+1,1,0,1)无法求解出正确答案,这个方程应该有两个解 i 和 -i,但变量数为1.
请点击"回复此楼",否则我将无法收到回帖提醒。
问题如果比较复杂或较难,请邮箱联系kimist@qq.com

论坛优秀回答者

权威

3630 麦片

财富积分



20

主题

3779

帖子

769

最佳答案
  • 关注者: 430
 楼主| 发表于 2017-3-20 16:21:06 | 显示全部楼层
本帖最后由 wuyou136 于 2017-3-20 16:24 编辑
kastin 发表于 2017-3-20 11:07
x,feval指的是第一个解吗?既然如此,为何不直接给multiSol和multiFeval呢?

刚刚试了复数解,不知道为 ...

x, fval 指的是 multiSol multiFval 中最好的解,当不需要给出多解时,x 和 fval 会比较方便。

对于复数解,我之前大意了,我采用的目标函数是本身带复数的,例如:fun = @(x)x^2+i ,这个目标函数能成功求得复数解。 对于 fun = @(x)x^2+1,如果初值不给复数的话,lsqnonlin 就只会在实数域进行搜索(MultiStart 里面用的搜索算法就是 lsqnonlin),所以优化出来就只能得到 x = 0。而 MultiStart 又给不了复数初值,所以这个问题好像还一下不太好解决......

我想到一个方法,就是再引入一个变量,比如对于目标函数:,我修改成:,然后同时优化
  1. fun = @(x)x.^2+1;
  2. fun_complex = @(x)fun(x(1)+i*x(2)); %修改目标函数

  3. obj_fun = fun_complex;
  4. start_num = 100;
  5. n = 1;
  6. lb = -inf; ub = inf;

  7. n = n*2; lb = repmat(lb,1,n); ub = repmat(ub,1,n);
  8. options = optimoptions('particleswarm','display','off');
  9. x0 = particleswarm(@(x)sum(obj_fun(x).^2),n,lb,ub,options);
  10. opts = optimoptions(@lsqnonlin ,'Display','off','Algorithm','Levenberg-Marquardt','FunValCheck','off');
  11. problem  = createOptimProblem('lsqnonlin','x0',x0,'objective',obj_fun,'lb',lb,'ub',ub,'options',opts);
  12. gs = MultiStart('FunctionTolerance',1e-1,'XTolerance',1e-1);
  13. [x,~,exitflag,~,allmins] = gs.run(problem,start_num);
  14. fval = obj_fun(x);
  15. x = x(1)+i*x(2) %计算解的实部和虚部
  16. fval

  17. function [c ceq] = constrain(x,f)
  18. ceq = f(abs(x));
  19. c =[];
  20. end
复制代码
  1. x =

  2.    0.0000 - 1.0000i


  3. fval =

  4.   -1.1386e-06 - 1.9404e-07i
复制代码

但是这样就多引入了变量,对于多元函数肯定大大降低搜索能力......所以对于复数解,还得想想有没有什么比较好的方法才行

论坛优秀回答者

退役版主

2206 麦片

财富积分



46

主题

2115

帖子

303

最佳答案
  • 关注者: 160
发表于 2017-3-20 20:06:08 | 显示全部楼层
wuyou136 发表于 2017-3-20 16:21
x, fval 指的是 multiSol multiFval 中最好的解,当不需要给出多解时,x 和 fval 会比较方便。

对于复数 ...

主要是因为复数无法比较大小,导致无法用区间lb,ub给出,因此直接求解可能无法行得通,不过似乎fsolve能自适应,应该可以想办法让MultiStart来获取复平面采样而不是数轴采样。
fsolve(@(x)x^2+1,1+i)

另外还有个问题就是,如果我们想要求非零解(很多问题中零解是平凡解),但lb和ub都是闭区间的两端,所以直接设置lb=0不行,稍微取小一点比如lb=eps,也是不行,因为数值求解的原因导致无论lb取什么值,只要不超过最小非零解(但问题在于根本不确定这个最小非零解的值是多少),返回的解都是非常接近lb的值(被零解强烈吸引)。不知道有什么好办法来解决?
请点击"回复此楼",否则我将无法收到回帖提醒。
问题如果比较复杂或较难,请邮箱联系kimist@qq.com

论坛优秀回答者

权威

3630 麦片

财富积分



20

主题

3779

帖子

769

最佳答案
  • 关注者: 430
 楼主| 发表于 2017-3-20 20:23:54 | 显示全部楼层
kastin 发表于 2017-3-20 20:06
主要是因为复数无法比较大小,导致无法用区间lb,ub给出,因此直接求解可能无法行得通,不过似乎fsolve能 ...

我觉得求复数解似乎应该不是很难,但是现在没有找到合适的方法,等看后面多读读 MultiStart 的文档看能不能找到解决方法。

你说的非零解的话,我的想法是,能不能利用 GlobalSolve 求得多解,然后在求得的解当中,排除掉第一个零解?
  1. f = @(x)sin(x)
  2. [x,fval,~,ms,mf] = GlobalSolve(f,1,eps,10)
复制代码

ms =

  1 至 2 列

   0.000000644724328   3.141592653589788

  3 至 4 列

   6.283185306824658   9.424777960769008
得到这4个解之后,再排除第一个解

论坛优秀回答者

退役版主

2206 麦片

财富积分



46

主题

2115

帖子

303

最佳答案
  • 关注者: 160
发表于 2017-3-20 21:28:36 | 显示全部楼层
wuyou136 发表于 2017-3-20 20:23
我觉得求复数解似乎应该不是很难,但是现在没有找到合适的方法,等看后面多读读 MultiStart 的文档看能不 ...

问题在于你这个方程很简单,多解就是固定的,对于复杂方程,会产生几十个解(但实际上真正的解只有4个)。算出来的多解中接近零的的解有很多个,但事实上他们都是因为初值的原因造成的假象,因为真正精确解中零解是一个平凡解(引力太大,将多个初始值都吸引下来了),剩下的非零解分散各处。现实中我们由于不知道精确的非零解与这些假的零解有多接近,造成不知道该排除哪些而不应该排除哪些。不知道我这样描述你是否理解?

请点击"回复此楼",否则我将无法收到回帖提醒。
问题如果比较复杂或较难,请邮箱联系kimist@qq.com

论坛优秀回答者

权威

3630 麦片

财富积分



20

主题

3779

帖子

769

最佳答案
  • 关注者: 430
 楼主| 发表于 2017-3-21 08:57:28 | 显示全部楼层
本帖最后由 wuyou136 于 2017-3-21 09:00 编辑
kastin 发表于 2017-3-20 21:28
问题在于你这个方程很简单,多解就是固定的,对于复杂方程,会产生几十个解(但实际上真正的解只有4个) ...

我明白你说的意思了,但是还是不太弄得懂造成这种情况的根本原因,如果我将非零解进行缩放:



通过控制 n 来使非零解和零解分离开来是否可行呢,你有没有这种问题的实际例子?

论坛优秀回答者

退役版主

2206 麦片

财富积分



46

主题

2115

帖子

303

最佳答案
  • 关注者: 160
发表于 2017-3-21 13:58:31 | 显示全部楼层
wuyou136 发表于 2017-3-21 08:57
我明白你说的意思了,但是还是不太弄得懂造成这种情况的根本原因,如果我将非零解进行缩放:

这种方法本质上是通过缩放x轴来拉开根之间的相对距离,但还是存在一些问题。

首先,因为方程中x缩小n倍,这意味着非零根之间也拉开了距离,然而求解时间与搜索空间是正相关的,因此这种做法牵一发而动全身,得不偿失。即使非零根之间距离非常近,只要与零点很接近也会难以区别。对于多元情形(比如二元f(x,y)=0,g(x,y)=0),这种方法更是难以凑效,特别是f的值随x,y变化快而g的值随x,y变化慢的情况,导致无法协调让二者同时到零的速度。

我有几个思路,不知是否总有效(以下均不考虑函数有间断点情形):
1. 如果方程不存在无穷远点的根,且非零根与零点距离大于1,可用倒代换消去零点影响;
2. 如果方程不存在无穷远点的根,但非零根附近f随某变量x变化振幅小,可f整体除以相应x;
3. 若同时有零点和无穷远点根,若能通过提出相应项(如x, 1/x, exp(-x)等)而消去这两个特殊根,则先手工处理,若不能明显消除则尝试乘以(x+1/x)这一项(多元情形类推)。

实际例子的表达式很难构造,但可以用文字模拟描述曲线变化情形。

以一元情况来说明,设曲线从零点穿过(可以是振荡穿过也可是平滑穿过)往上升,上升到波峰(波峰高度大概在1e-3量级)后下降,下降到1e-5高度开始高频小幅度振荡(振幅量级大概1e-5~1e-6),右行然后平缓逼近x轴,与x轴相切(非零点),然后平缓右行,再次小幅度高频振荡(如前,振荡曲线始终可以保持少数次穿过x轴,也可不穿过),平缓穿过x轴后还可继续振荡或者始终近似平行于x轴延伸……

这一切都是发生在[0,1]区间内,且非常接近零点(比如发生在[0,1e-5]),这种形态的函数曲线存在非常多个局部极值(非零根之间也相距紧密),导致全局搜索容易陷入局部最优。所以无法分辨到底求出来的多解是真实的根还只是因为精度和算法造成的局部最优解(假根)。

上面的例子太过于复杂,很难找到合适的表达式满足要求,不过这里可以构造一个简单的例子(虽然表达式上容易直接剥离x^2项,但这里我们只是等价模拟,即有其他复杂的表达式能跟该函数曲线极其接近,但无从从表达式上看出这一点)试试,不知能否说明:
f=@(x)0.7886525862862048*x^10-0.9487490613023043*x^9-0.365551290385882*x^8+0.4431736989549043*x^7+0.06077675574420972*x^6-0.07470986247652193*x^5+0.0002235706607776218*x^4+1.371386992637875*10^-25*x^3-4.103901224851498*10^-28*x^2
方程的根精确根为 x=-1.35484946596096*10^-12, x=0(二重根),x=1.35484946596096*10^-12, x=0.003, x=1.2
请点击"回复此楼",否则我将无法收到回帖提醒。
问题如果比较复杂或较难,请邮箱联系kimist@qq.com

论坛优秀回答者

权威

3630 麦片

财富积分



20

主题

3779

帖子

769

最佳答案
  • 关注者: 430
 楼主| 发表于 2017-3-22 16:04:07 | 显示全部楼层
kastin 发表于 2017-3-21 13:58
这种方法本质上是通过缩放x轴来拉开根之间的相对距离,但还是存在一些问题。

首先,因为方程中x缩小n倍 ...

谢谢,解释的很详细,完全明白其中的难点所在了。我试了下坐标轴缩放方法,对于你提供的例子,手动不断调节x轴和y轴的缩放才能勉强把第一个非零根找到,但是用这种方法自动把所有根找齐不太可能,对于多元函数就更困难了。

我数学功底太差了,现在没有什么好的想法,后面我会一直关注这个问题,如果有好的想法马上找你商量。

论坛优秀回答者

退役版主

2206 麦片

财富积分



46

主题

2115

帖子

303

最佳答案
  • 关注者: 160
发表于 2017-3-22 18:54:39 | 显示全部楼层
wuyou136 发表于 2017-3-22 16:04
谢谢,解释的很详细,完全明白其中的难点所在了。我试了下坐标轴缩放方法,对于你提供的例子,手动不断调 ...

好的,这是因为前些天想用matlab求解一个多元高次多项式方程组的所有解,想到看用你的的这个函数能否行得通,结果恰好就是产生这种问题。从全局优化来看,这个问题本身还是跟初值有大大的关系,初值如果选择的好则直接收敛到真解;如何消除初始值带来的问题,看来只能用遗传算法之类的了,因为它有可能跳出局部最优。
请点击"回复此楼",否则我将无法收到回帖提醒。
问题如果比较复杂或较难,请邮箱联系kimist@qq.com
您需要登录后才可以回帖 登录 | 注册

本版积分规则

关闭

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

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