楼主: wuyou136

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

  [复制链接]

论坛优秀回答者

中级

1040 麦片

财富积分


5001500


0

主题

2545

帖子

225

最佳答案
  • 关注者: 146
发表于 2017-3-24 10:48:48 | 显示全部楼层
手头没有安装Matlab,无法享用楼主的方程组求解工具箱,哪位试试解下面的方程组(仅限实数解),结果如何?

6250*pi*x1+2500*pi*x1*(15-x2)+211/20-x3=0;
-169169.3148*x2+168971.8978*x2-x4-1221709=0;
1.03*x2+1.84*x2*x4+2500*x4-x3=0;
6166.837109-282.1863256*x2-x1*(20-x2)=0;

论坛优秀回答者

权威

3630 麦片

财富积分



20

主题

3779

帖子

769

最佳答案
  • 关注者: 431
 楼主| 发表于 2017-3-24 14:51:57 | 显示全部楼层
shihe 发表于 2017-3-24 10:48
手头没有安装Matlab,无法享用楼主的方程组求解工具箱,哪位试试解下面的方程组(仅限实数解),结果如何? ...

这个方程组按默认参数的话得不到收敛解,要将搜索初值数给的很大才能得到收敛解,这里我给到了3000,不过同时计算时间也非常长,我开了并行计算还计算了20min才得到结果:
  1. syms x1 x2 x3 x4
  2. eq1 = 6250*pi*x1+2500*pi*x1*(15-x2)+211/20-x3 ;
  3. eq2 = -169169.3148*x2+168971.8978*x2-x4-1221709 ;
  4. eq3 = 1.03*x2+1.84*x2*x4+2500*x4-x3 ;
  5. eq4 = 6166.837109-282.1863256*x2-x1*(20-x2) ;
  6. eq = matlabFunction([eq1,eq2,eq3,eq4],'vars',{[x1,x2,x3,x4]})

  7. tic
  8. [x,fval,exitflag,ms,mf] = GlobalSolve(eq,4,[],[],3e3);
  9. toc
  10. x.'
  11. fval.'
复制代码
  1. x =
  2.    1.0e+09 *

  3.    0.000158561829777
  4.    0.000000019996695
  5.   -3.109238410946461
  6.   -0.001225656687541

  7. fval =
  8. 1.0e-06 *

  9.   -0.190734862570707
  10.                    0
  11.    0.111758708953857
  12.   -0.000060936145019
复制代码

主要原因是收敛解达到了10^9数量级,搜索空间越大,耗费的时间也越大。
另外,这个方程其实是有解析解的,直接用 vpasolve 也能得到结果:
  1. syms x1 x2 x3 x4 real
  2. eq1 = 6250*pi*x1+2500*pi*x1*(15-x2)+211/20-x3 ;
  3. eq2 = -169169.3148*x2+168971.8978*x2-x4-1221709 ;
  4. eq3 = 1.03*x2+1.84*x2*x4+2500*x4-x3 ;
  5. eq4 = 6166.837109-282.1863256*x2-x1*(20-x2) ;
  6. [x1 x2 x3 x4] = vpasolve([eq1,eq2,eq3,eq4],[x1 x2 x3 x4])
复制代码
  1. x1 =

  2. 158561.82977690677275591023419543


  3. x2 =

  4. 19.996695022900017268780308373052


  5. x3 =

  6. -3109238410.9464602771027430761724


  7. x4 =

  8. -1225656.6875413355873260785336174
复制代码

新手

5 麦片

财富积分


050


3

主题

27

帖子

0

最佳答案
发表于 2017-5-16 20:23:42 | 显示全部楼层
wuyou136 发表于 2017-2-26 20:11
需要直接使用的看二楼:
将 GlobalSolve 下载放到搜索路径之后,输入:
修改你的目标函数即可
  1. 错误使用 optimoptions (line 114)
  2. Invalid solver specified. Provide a solver name or handle (such as 'fmincon' or @fminunc).
  3. Type DOC OPTIMOPTIONS for a list of solvers.

  4. 出错 GlobalSolve (line 28)
  5. options = optimoptions('particleswarm','display','off');
复制代码

为什么我的是这样啊,我已经将GlobalSolve添加到路径下了啊

论坛优秀回答者

权威

3630 麦片

财富积分



20

主题

3779

帖子

769

最佳答案
  • 关注者: 431
 楼主| 发表于 2017-5-16 20:47:55 | 显示全部楼层
周周迟到 发表于 2017-5-16 20:23
为什么我的是这样啊,我已经将GlobalSolve添加到路径下了啊

应该是版本问题,建议升级 MATLAB 版本

新手

5 麦片

财富积分


050


3

主题

27

帖子

0

最佳答案
发表于 2017-5-18 19:25:34 | 显示全部楼层
wuyou136 发表于 2017-5-16 20:47
应该是版本问题,建议升级 MATLAB 版本
  1. clear;
  2. clc;

  3. % % 参数输入
  4. t=0.02;%板厚
  5. r0=2.5;%初始曲率半径
  6. r=1.9;%目标半径

  7. % kfi_1=1/(r0+0.5*t);%初始曲率kfi-1
  8. kfi_1=0
  9. kfi=1/(r+0.5*t);%目标的曲率

  10. u=0;%摩擦系数
  11. E=2.06*10^11;%平面应变的杨氏模量
  12. I=1/12*t.^3;
  13. E1=2726*10^6;%切线模量
  14. re=E1/E;%硬化程度E1/E
  15. Y=235*10^6;%屈服极限
  16. c=(3.^0.5*E*(1./(1-0.28.^2)-4./3*re))./(2*(1-re)*Y);
  17. kei=2./(t*c)+(c+1)./c*kfi_1;
  18. ke1=2./(t*c);

  19. % 计算k2i
  20. syms k y;
  21. yep=1./(c*k-(c+1)*kfi_1);%屈服点到中性层的距离
  22. Y1=E./(1-0.28.^2)*(k-kfi_1)*y./(1+kfi_1*y);%弹性阶段的应力表达式
  23. Y2=2*3.^(-0.5)*(1-re)*Y+4./3*E1*(k-kfi_1)*y./(1+kfi_1*y);%塑性阶段的应力表达式
  24. M1=2*int(Y1*y,y,0,0.5*t);%弹性阶段弯矩与曲率关系
  25. M2=2*int(Y1*y,y,0,yep)+2*int(Y2*y,y,0,0.5*t)-2*int(Y2*y,y,0,yep);%塑性阶段弯矩与曲率关系
  26. fun=k-(1-0.28.^2)*M2./(E*I)-kfi;
  27. fun=matlabFunction(fun);
  28. x0=7;
  29. options=optimset('TolFun',1e-11,'TolX',1e-11);
  30. [n,fval,exitflag,output]=fzero(fun,x0,options);
  31. k2i=n;

  32. % % 现在开始计算O1、O2
  33. % 卷板机基本参数
  34. r1=0.475;%下辊半径0.475
  35. r2=0.54;%上辊半径
  36. r3=0.4;%侧辊半径0.4
  37. R1=r1+0.5*t;
  38. R2=r2+0.5*t;
  39. R3=r3+0.5*t;
  40. Kp=0.5*pi*212*10^9./(1-0.28^2);
  41. w=atan(u);%摩擦角
  42. D=tan(25/180*pi);
  43. x01=(1.27577+r2+r3+t)*D;%侧辊初始横坐标
  44. m=(r2+0.5*t)./(r1+0.5*t);%中性层半径比
  45. B=0.5*E*I*(k2i.^2-kfi.^2)./(1-0.28.^2);
  46. M22=E*I*(k2i-kfi)./(1-0.28.^2);

  47. syms O2 O3 k;
  48. T1=O3-asin((k.^2-kfi.^2)*sin(O3-O2)./(k2i.^2-kfi.^2));%x=O3 y=k,f(k)的表达式
  49. t1=diff(T1,k,1);
  50. f1=cos(T1)*t1./k;
  51. f2=sin(T1)*t1./k;
  52. f1=matlabFunction(f1);
  53. f2=matlabFunction(f2);
  54. h1=@(O2,O3)(integral(@(k)f1(O2,O3,k),k2i,kfi))*cos(O3);%(x3-x2)*cos(O3)
  55. h2=@(O2,O3)(integral(@(k)f2(O2,O3,k),k2i,kfi))*sin(O3);%(y3-y2)*sin(O3)
  56. eq1=@(O2,O3)(0.5*E*I*(k2i.^2-kfi.^2)./(1-0.28.^2)./sin(O3-O2))*(h1(O2,O3)+h2(O2,O3))-M22;
  57. eq2=@(O2,O3)(r3+0.5*t)*(sin(O3)+D*(1-cos(O3)))+(r2+0.5*t)*(sin(O2)+D*(1-cos(O2)))+integral(@(k)(f1(O2,O3,k)+D*f2(O2,O3,k)),k2i,kfi)-x01;
  58. % options=optimset('Algorithm','levenberg-marquardt','Display','iter','TolFun',1e-11,'TolX',1e-11,'FunValCheck','on');
  59. eq=@(x)[eq1(x(1),x(2));eq2(x(1),x(2))];
  60. [x,fval,exitFlag,multiSol] = GlobalSolve1(eq,2)
复制代码

这是我的问题计算过程,要求O2、O3,根据您提供的代码一下子求出了99个解,但是大部分值是我不需要的,我大概知道O2的范围(0,0.03),O3的范围(0,0.7),请问我怎么设置搜索区间啊,谢谢。

论坛优秀回答者

权威

3630 麦片

财富积分



20

主题

3779

帖子

769

最佳答案
  • 关注者: 431
 楼主| 发表于 2017-5-18 19:30:15 | 显示全部楼层
周周迟到 发表于 2017-5-18 19:25
这是我的问题计算过程,要求O2、O3,根据您提供的代码一下子求出了99个解,但是大部分值是我不需要的, ...
  1. [x,fval,exitFlag,multiSol] = GlobalSolve1(eq,2,[0, 0], [0.3, 0.7])
复制代码

x =

    0.0819    0.5200

新手

5 麦片

财富积分


050


3

主题

27

帖子

0

最佳答案
发表于 2017-5-18 19:35:11 | 显示全部楼层
wuyou136 发表于 2017-5-18 19:30
x =

    0.0819    0.5200

十分感谢,请问一下返回值x只是其中的一个解multiSol{1}吗?还是x有其他含义

论坛优秀回答者

权威

3630 麦片

财富积分



20

主题

3779

帖子

769

最佳答案
  • 关注者: 431
 楼主| 发表于 2017-5-18 19:36:09 | 显示全部楼层
周周迟到 发表于 2017-5-18 19:35
十分感谢,请问一下返回值x只是其中的一个解multiSol{1}吗?还是x有其他含义 ...

x 代表的是 multiSol 中最好的解

新手

5 麦片

财富积分


050


3

主题

27

帖子

0

最佳答案
发表于 2017-5-18 19:48:52 | 显示全部楼层
wuyou136 发表于 2017-5-18 19:36
x 代表的是 multiSol 中最好的解

我发现什么参数都不改,每次计算结果有点不一样,请问为什么啊

论坛优秀回答者

权威

3630 麦片

财富积分



20

主题

3779

帖子

769

最佳答案
  • 关注者: 431
 楼主| 发表于 2017-5-18 19:51:29 | 显示全部楼层
周周迟到 发表于 2017-5-18 19:48
我发现什么参数都不改,每次计算结果有点不一样,请问为什么啊

因为这里采用的是随机初值进行搜索,如果方程有多解的话,就有可能每次收敛的解不一样
您需要登录后才可以回帖 登录 | 注册

本版积分规则

关闭

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

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