楼主: wuyou136

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

  [复制链接]

新手

5 麦片

财富积分


050


3

主题

27

帖子

0

最佳答案
发表于 2017-5-19 10:24:09 | 显示全部楼层
wuyou136 发表于 2017-5-18 19:51
因为这里采用的是随机初值进行搜索,如果方程有多解的话,就有可能每次收敛的解不一样 ...

如果是一个一元方程,要求所有的解,该如何在这个基础上修改啊

论坛优秀回答者

权威

3630 麦片

财富积分



20

主题

3779

帖子

769

最佳答案
  • 关注者: 431
 楼主| 发表于 2017-5-19 10:25:41 | 显示全部楼层
周周迟到 发表于 2017-5-19 10:24
如果是一个一元方程,要求所有的解,该如何在这个基础上修改啊

直接套用函数就行,ms 就是求得的多解

新手

5 麦片

财富积分


050


3

主题

27

帖子

0

最佳答案
发表于 2017-5-19 10:30:23 | 显示全部楼层
wuyou136 发表于 2017-5-19 10:25
直接套用函数就行,ms 就是求得的多解
  1. clear;
  2. clc;

  3. % % 参数输入
  4. t=0.02;%板厚
  5. r0=3.4;%初始曲率半径
  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. D=tan(25/180*pi);
  41. x01=(1.27577+r2+r3+t)*D;%侧辊初始横坐标
  42. m=(r2+0.5*t)./(r1+0.5*t);%中性层半径比
  43. B=0.5*E*I*(k2i.^2-kfi.^2)./(1-0.28.^2);
  44. M22=E*I*(k2i-kfi)./(1-0.28.^2);

  45. syms O2 O3 k;
  46. T1=O3-asin((k.^2-kfi.^2)*sin(O3-O2)./(k2i.^2-kfi.^2));%x=O3 y=k,f(k)的表达式
  47. t1=diff(T1,k,1);
  48. f1=cos(T1)*t1./k;
  49. f2=sin(T1)*t1./k;
  50. f1=matlabFunction(f1);
  51. f2=matlabFunction(f2);
  52. h1=@(O2,O3)(integral(@(k)f1(O2,O3,k),k2i,kfi))*cos(O3);%(x3-x2)*cos(O3)
  53. h2=@(O2,O3)(integral(@(k)f2(O2,O3,k),k2i,kfi))*sin(O3);%(y3-y2)*sin(O3)
  54. 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;
  55. 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;
  56. % options=optimset('Algorithm','levenberg-marquardt','Display','iter','TolFun',1e-11,'TolX',1e-11,'FunValCheck','on');
  57. eq=@(x)[eq1(x(1),x(2));eq2(x(1),x(2))];
  58. [x,fval,exitFlag,multiSol] = GlobalSolve1(eq,2,[0,0.2], [0.1,0.8])  %用hua和myself的方法发现成形半径越小,O3越大,O3的范围(0.3,0.8)O2很小,几乎不超过0.1
  59. % % 之前O2、O3的范围[0,0.3], [0.1,0.8]
  60. max=size(multiSol,2)%查看数组multiSol有多少列

  61. % % 取O2、O3的最优值
  62. for i=1:max;  
  63. % [A1 A2 A3 A4]=multiSol{1:4}
  64. O2=multiSol{i}(1)
  65. theat2(i)=O2
  66. O3=multiSol{i}(2)
  67. theat3(i)=O3
  68. xs=(r3+0.5*t)*sin(O3)+(r2+0.5*t)*sin(O2)+integral(@(k)f1(O2,O3,k),k2i,kfi);
  69. ys=-(r3+0.5*t)*cos(O3)+(r2+0.5*t)*(1-cos(O2))+integral(@(k)f2(O2,O3,k),k2i,kfi);
  70. format short;
  71. hx=(x01-xs)*1000  %水平位移增量
  72. HX(i)=hx  %将hx的值存储到向量HX中
  73. vy=(r3+0.5*t+ys)*1000  %竖直位移增量
  74. VY(i)=vy
  75. s=(hx.^2+vy.^2).^0.5  %总的进给量
  76. S(i)=s
  77. end


  78. O2=0.0118
  79. % 计算力的大小
  80. fun2=@(x)(B./sin(O2-x))*((R2*sin(O2)+R1*sin(x))*cos(x)+(R2*(1-cos(O2))+R1*(1-cos(x)))*sin(x))-M22;
  81. % options=optimset('TolFun',1e-11,'TolX',1e-11);
  82. % [n,fval,exitflag,output]=fzero(fun2,0.005,options);
  83. [x,fval,exitFlag,multiSol] = GlobalSolve1(fun2,1,[0,0.2])
复制代码


O2是我根据您的公式求出的解,然后我想在用全局搜索找fun2的解,但是运行之后会报错
未定义函数或变量 'lb'。

出错 GlobalSolve1 (line 29)
x0 = particleswarm(@(x)sum(obj_fun(x).^2),n,lb,ub,options);
请问这是怎么回事

论坛优秀回答者

权威

3630 麦片

财富积分



20

主题

3779

帖子

769

最佳答案
  • 关注者: 431
 楼主| 发表于 2017-5-19 10:33:25 | 显示全部楼层
周周迟到 发表于 2017-5-19 10:30
O2是我根据您的公式求出的解,然后我想在用全局搜索找fun2的解,但是运行之后会报错
未定义函数或变量  ...

[x,fval,exitFlag,multiSol] = GlobalSolve1(fun2,1,0,0.2)
上下限要分开

新手

5 麦片

财富积分


050


0

主题

12

帖子

0

最佳答案
发表于 2017-5-23 11:35:13 | 显示全部楼层
楼主您好!我的方程组用这种方法求解一致出现问题,不知道如何解决,请求楼主帮我测试下,我的方程组如下:
function F=myfun(x)
x0=x(1);
y0=x(2);
z0=x(3);
A=x(4);
B=x(5);
D=pi/6;
x1=0 ; y1=0;  z1=0;
x2=0 ; y2=4;  z2=0;
x3=2.5 ; y3=0;  z3=0;
x4=5 ; y4=0;  z4=0;
x5=5 ; y5=2;  z5=0;
x6=0 ; y6=5;  z6=0;
t1=0.006784835;
t2=0.00390225;
t3=0.006102243;
t4=0.009790904;
t5=0.00723989;
t6=0.005564238;
F(1)=((x1-(680*t1*cos(A)*sin(B)+x0))^2+(y1-(680*t1*sin(A)+y0))^2+(z1-(-680*t1*cos(A)*cos(B)+z0))^2)^0.5*cos(D)-(-(x1-(680*t1*cos(A)*sin(B)+x0))*cos(A)*sin(B)-(y1-(680*t1*sin(A)+y0))*sin(A)+(z1-(-680*t1*cos(A)*cos(B)+z0))*cos(A)*cos(B));
F(2)=((x2-(680*t2*cos(A)*sin(B)+x0))^2+(y2-(680*t2*sin(A)+y0))^2+(z2-(-680*t2*cos(A)*cos(B)+z0))^2)^0.5*cos(D)-(-(x2-(680*t2*cos(A)*sin(B)+x0))*cos(A)*sin(B)-(y2-(680*t2*sin(A)+y0))*sin(A)+(z2-(-680*t2*cos(A)*cos(B)+z0))*cos(A)*cos(B));
F(3)=((x3-(680*t3*cos(A)*sin(B)+x0))^2+(y3-(680*t3*sin(A)+y0))^2+(z3-(-680*t3*cos(A)*cos(B)+z0))^2)^0.5*cos(D)-(-(x3-(680*t3*cos(A)*sin(B)+x0))*cos(A)*sin(B)-(y3-(680*t3*sin(A)+y0))*sin(A)+(z3-(-680*t3*cos(A)*cos(B)+z0))*cos(A)*cos(B));
F(4)=((x4-(680*t4*cos(A)*sin(B)+x0))^2+(y4-(680*t4*sin(A)+y0))^2+(z4-(-680*t4*cos(A)*cos(B)+z0))^2)^0.5*cos(D)-(-(x4-(680*t4*cos(A)*sin(B)+x0))*cos(A)*sin(B)-(y4-(680*t4*sin(A)+y0))*sin(A)+(z4-(-680*t4*cos(A)*cos(B)+z0))*cos(A)*cos(B));
F(5)=((x5-(680*t5*cos(A)*sin(B)+x0))^2+(y5-(680*t5*sin(A)+y0))^2+(z5-(-680*t5*cos(A)*cos(B)+z0))^2)^0.5*cos(D)-(-(x5-(680*t5*cos(A)*sin(B)+x0))*cos(A)*sin(B)-(y5-(680*t5*sin(A)+y0))*sin(A)+(z5-(-680*t5*cos(A)*cos(B)+z0))*cos(A)*cos(B));
F(6)=((x6-(680*t6*cos(A)*sin(B)+x0))^2+(y6-(680*t6*sin(A)+y0))^2+(z6-(-680*t6*cos(A)*cos(B)+z0))^2)^0.6*cos(D)-(-(x6-(680*t6*cos(A)*sin(B)+x0))*cos(A)*sin(B)-(y6-(680*t6*sin(A)+y0))*sin(A)+(z6-(-680*t6*cos(A)*cos(B)+z0))*cos(A)*cos(B));

论坛优秀回答者

权威

3630 麦片

财富积分



20

主题

3779

帖子

769

最佳答案
  • 关注者: 431
 楼主| 发表于 2017-5-23 16:17:13 | 显示全部楼层
li329895720 发表于 2017-5-23 11:35
楼主您好!我的方程组用这种方法求解一致出现问题,不知道如何解决,请求楼主帮我测试下,我的方程组如下: ...

5个未知数,6个方程一般都是无解的吧,我用 GlobalSolve 没有找到收敛解

新手

5 麦片

财富积分


050


0

主题

12

帖子

0

最佳答案
发表于 2017-5-23 16:55:54 | 显示全部楼层
有一个方程式是多的,任意5个方程组可以应该解出一组相同的解,楼主可以把使用GlobalSolve函数的文件发我下吗?我调用一直不成功。

新手

5 麦片

财富积分


050


0

主题

12

帖子

0

最佳答案
发表于 2017-5-23 16:56:13 | 显示全部楼层
wuyou136 发表于 2017-5-23 16:17
5个未知数,6个方程一般都是无解的吧,我用 GlobalSolve 没有找到收敛解


有一个方程式是多的,任意5个方程组可以应该解出一组相同的解,楼主可以把使用GlobalSolve函数的文件发我下吗?我调用一直不成功

论坛优秀回答者

权威

3630 麦片

财富积分



20

主题

3779

帖子

769

最佳答案
  • 关注者: 431
 楼主| 发表于 2017-5-23 17:08:34 | 显示全部楼层
li329895720 发表于 2017-5-23 16:56
有一个方程式是多的,任意5个方程组可以应该解出一组相同的解,楼主可以把使用GlobalSolve函数的文件发 ...

我用的就是7L的文件

新手

5 麦片

财富积分


050


0

主题

12

帖子

0

最佳答案
发表于 2017-5-23 17:39:36 | 显示全部楼层
wuyou136 发表于 2017-5-23 17:08
我用的就是7L的文件

我调用的也是7楼文件,方程形式改为下式,还是提醒有错误!
syms t1  x0 y0 z0  A B D x y z;
eq1=((x-(680*t1*cos(A)*sin(B)+x0))^2+(y-(680*t1*sin(A)+y0))^2+(z-(-680*t1*cos(A)*cos(B)+z0))^2)^0.5*cos(D)-(-(x-(680*t1*cos(A)*sin(B)+x0))*cos(A)*sin(B)-(y-(680*t1*sin(A)+y0))*sin(A)+(z-(-680*t1*cos(A)*cos(B)+z0))*cos(A)*cos(B));
eq11= subs(eq1,[D x y z t1], [pi/6 0 0 0 0.001971118]);       
eq21= subs(eq1,[D x y z t1], [pi/6 0 5 0 0.007317894]);       
eq31= subs(eq1,[D x y z t1], [pi/6 2.5 0 0 0.003680043]);               
eq41= subs(eq1,[D x y z t1], [pi/6 5 0 0 0.010301765]);
eq51= subs(eq1,[D x y z t1], [pi/6 5 2.5 0 0.010318838]);
eq = matlabFunction([eq11,eq21,eq31,eq41,eq51],'vars',{[x0,y0,z0,A,B]})
tic
[x,fval,exitflag,ms,mf] = GlobalSolve(eq,4,[],[],3e3);
toc
x.'
fval.'
您需要登录后才可以回帖 登录 | 注册

本版积分规则

关闭

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

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