全局最优化的讨论

rocwoods

我们知道，1stOpt具有全局最优的能力，可惜其算法无缘一见。全局最优是个多年的难题，其方法目前有好多种。开本帖的目的就是希望各路高手能参与讨论，我们虽然未必一定能开创新的有效算法，但是如果我们能够将现有的多种算法结合算例展现在这里，那也是一个全局最优化不错的总结贴，也不枉讨论！
      言归正传，第一个开刀的例子就是1stOpt帮助文档里给的针状函数，求其全局最大值：

rocwoods

下面给出我以前编写的一个随机行走法求函数最优值的方法，函数基本代码如下：

1. 
   
2. function [mx,minf]=randwalk(f,x,lamda,epsilon,N)
   
3. %随机行走法求函数的极小值。输入f为所求函数的句柄，
   
4. %x为初始值。lamda为步长。epsilon为控制lamda的减小的阈值，即lamda减小到epsilon时迭代停止。
   
5. %N为为了产生较好点的迭代控制次数。
   
6. %函数返回值mx为n次试验得到的最优解，minf为响应的最优值。
   
7. %copyright: rocwoods
   
8. f1=f(x(1),x(2));
   
9. while(lamda>=epsilon)
   
10.     k=1;
    
11.     while(k<=N)
    
12.         u=unifrnd(-1,1,1,2);
    
13.         x1=x+lamda*(u/norm(u));
    
14.         f11=f(x1(1),x1(2));
    
15.         if f11<f1
    
16.             f1=f11;
    
17.             x=x1;
    
18.             k=1;
    
19.         else
    
20.             k=k+1;
    
21.         end
    
22.     end
    
23.     lamda=lamda/2;
    
24. end
    
25. mx=x1;
    
26. minf=f1;
    

复制代码

事实证明，当[x,y]位于以(50,50)为圆心，半径为1点多的圆域内，randwalk法几乎都能找到最优点。下面给出该方法和MATLAB自带的fminsearch的对比:
比较前的准备工作，建立各自需要的目标函数的函数句柄：

1. 
   
2. f=@(x) -sin(sqrt( (x(1)-50)^2+(x(2)-50)^2 ) +exp(1))/(sqrt( (x(1)-50)^2+(x(2)-50)^2 ) +exp(1))-1;
   
3. f2=@(x,y) -sin(sqrt( (x-50).^2+(y-50).^2 ) +exp(1))./(sqrt( (x-50).^2+(y-50).^2 ) +exp(1))-1;
   

复制代码

比较：

1. 
   
2. >> [x,fval] = fminsearch(f, [49,49], optimset('TolX',1e-8))
   
3. x =
   
4.    51.4416   45.2051
   
5. fval =
   
6.    -1.1284
   
7. >> [x,fval] = fminsearch(f, [48.5,48.5], optimset('TolX',1e-8))
   
8. x =
   
9.    48.1592   45.3437
   
10. fval =
    
11.    -1.1284
    
12. >> [x,fval] = fminsearch(f, [48.9,48.9], optimset('TolX',1e-8))
    
13. x =
    
14.    50.8189   45.0605
    
15. fval =
    
16.    -1.1284
    
17. >> [x,fval] = fminsearch(f, [49.3,49.3], optimset('TolX',1e-8))
    
18. x =
    
19.    50.4007   45.0091
    
20. fval =
    
21.    -1.1284
    
22. >> [x,fval] = fminsearch(f, [49.5,49.5], optimset('TolX',1e-8))
    
23. x =
    
24.    51.3046   45.1660
    
25. fval =
    
26.    -1.1284
    
27. >> [x,fval] = fminsearch(f, [49.6,49.6], optimset('TolX',1e-8))
    
28. x =
    
29.    51.6997   45.2904
    
30. fval =
    
31.    -1.1284
    
32. >> [x,fval] = fminsearch(f, [49.7,49.7], optimset('TolX',1e-8))
    
33. x =
    
34.    50.0000   50.0000
    
35. fval =
    
36.    -1.1511
    
37. >> [x,fval] = fminsearch(f, [49.65,49.65], optimset('TolX',1e-8))
    
38. x =
    
39.    50.8348   45.0631
    
40. fval =
    
41.    -1.1284
    
42. >> [mx,minf]=randwalk(f2,[49,49],0.5,0.00001,100)
    
43. mx =
    
44.    50.0000   50.0000
    
45. minf =
    
46.    -1.1511
    
47. >> [mx,minf]=randwalk(f2,[48.5,48.5],0.5,0.00001,100)
    
48. mx =
    
49.    46.3953   46.5250
    
50. minf =
    
51.    -1.1284
    
52. >> [mx,minf]=randwalk(f2,[48.9,48.9],0.5,0.00001,100)
    
53. mx =
    
54.    50.0000   50.0000
    
55. minf =
    
56.    -1.1511
    
57. >> [mx,minf]=randwalk(f2,[49.3,49.3],0.5,0.00001,100)
    
58. mx =
    
59.    50.0000   50.0000
    
60. minf =
    
61.    -1.1511
    

复制代码

可见，randwalk比fminsearch的有效搜索范围要大，fminsearch在初始值(49.65,49.65)时按如上设置都不能迭代到全局最优解(调整函数设置不知能否实现)。而随机行走法在初始值为(48.9,48.9)时也能迭代到最优。

[ 本帖最后由 rocwoods 于 2008-12-5 17:41 编辑 ]

FreddyMusic

I played one D random walk in Mathematica, 2D not yet.

shamohu

看看这个SPSA优化算法：http://www.jhuapl.edu/SPSA/

不知与随机行走法有无共同点。

rocwoods

看了一下SPSA算法，似乎还是跟函数梯度有关。这点跟随机行走法还是不一样。
     基础的随机行走法，其求解过程是这样的：首先限定初始迭代点x，初次行走步长lamda($$$)，以及为了产生较好点的迭代控制次数N。k置初值1，（***），while k<N （&&&）,随机生成一个(-1,1)之间的二维向量u,标准化u：u/norm(u)，并令x1=x+lamda*(u/norm(u));完成第一步行走，如果f(x1)小于f(x)，即找到一个比初始值好的点，那么k重新置1，回到(***)处。否则k=k+1。回到（&&&）处,如果连续N次都找不到更优的值，我们认为，最优解就在以当前最优解为中心，当前步长为半径的球内，这时候，我们减半lamda,回到($$$)处，开始新一轮行走。这个过程直到lamda<epsilon时候结束。
      在此基础上，我提出了自己的一些改进方法。即在（&&&）处，随机生成n个(-1,1)之间的二维向量u，而不是一个，这样对应n个x1。从这n个x1中选取使得f(x1)最小的x1作为行走目标。这样大大加大了全局寻优的能力！！而且对初始值的依赖大大降低。以下是改进代码以及寻优效果：

1. function [mx,minf]=randwalk(f,x,lamda,epsilon,N,n)
   
2. %随机行走法求函数的极小值。输入f为所求函数的句柄，
   
3. %x为初始值。lamda为步长。epsilon为控制lamda的减小的阈值，即lamda减小到epsilon时迭代停止。
   
4. %N为为了产生较好点的迭代控制次数。
   
5. %n为单步循环行走次数，目的是尽可能走到全局最优点附近
   
6. %函数返回值mx为n次试验得到的最优解，minf为响应的最优值。
   
7. F=zeros(n,1);
   
8. X=zeros(n,2);
   
9. f1=f(x(1),x(2));
   
10. while(lamda>=epsilon)
    
11.     k=1;
    
12.     while(k<=N)
    
13.         u=unifrnd(-1,1,n,2);
    
14.         for ii=1:n
    
15.         X(ii,:)=x+lamda*(u(ii,:)/norm(u(ii,:)));
    
16.         F(ii)=f(X(ii,1),X(ii,2));
    
17.         end
    
18.         [f11,kk]=min(F);
    
19.         if f11<f1
    
20.             f1=f11;
    
21.             x=X(kk,:);
    
22.             k=1;
    
23.         else
    
24.             k=k+1;
    
25.         end
    
26.     end
    
27.     lamda=lamda/2;
    
28. end
    
29. mx=X(kk,:);
    
30. minf=f1;

复制代码

效果：

1. [mx,minf]=randwalk(f2,[0,0],10,0.00001,1000,10)
   
2. mx =
   
3.    50.0000   50.0000
   
4. minf =
   
5.    -1.1511
   
6. >> [mx,minf]=randwalk(f2,[-10,6],10,0.00001,1000,10)
   
7. mx =
   
8.    50.0000   50.0000
   
9. minf =
   
10.    -1.1511
    
11. >> [mx,minf]=randwalk(f2,[2.9765,0.888],10,0.00001,1000,10)
    
12. mx =
    
13.    50.0000   50.0000
    
14. minf =
    
15.    -1.1511
    
16. >> [mx,minf]=randwalk(f2,[12.231,3.777],10,0.00001,1000,10)
    
17. mx =
    
18.    50.0000   50.0000
    
19. minf =
    
20.    -1.1511
    
21. >> [mx,minf]=randwalk(f2,[pi,exp(0.653)],10,0.00001,1000,10)
    
22. mx =
    
23.    50.0000   50.0000
    
24. minf =
    
25.    -1.1511

复制代码

[ 本帖最后由 rocwoods 于 2007-8-7 11:03 编辑 ]

rocwoods

        首先建立隐函数的句柄。为了方便使用者，代码中所有初值都是随机初值。

1. z=@(x,y) fzero(@(z) z-sin((z*x-0.5)^2 + x*2*y^2-z/10)*exp(-((x-0.5-exp(-y+z))^2 + y^2-z/5+3)),rand);

复制代码

大家从下面测试中可以看出改进的随机行走法的适用范围非常具有弹性。

1. tic,[mx,minf]=randwalk(z,[unifrnd(-1,7),unifrnd(-2,2)],5,0.000001,100,5),toc
   
2. mx =
   
3.    2.898269787201496  -0.857313567083922
   
4. minf =
   
5.   -0.023354083261443
   
6. Elapsed time is 18.690406 seconds.（由于x,y被限定在[-1，7]，[-2，2]，因此初值步长设为5左右和区间半径相仿，初值随机给出，每步有5个位置供挑选）
   
7. >> tic,[mx,minf]=randwalk(z,[unifrnd(-1,7),unifrnd(-2,2)],2.5,0.000001,100,5),toc
   
8. mx =
   
9.    2.898271635033464  -0.857313487189683
   
10. minf =
    
11.   -0.023354083261454
    
12. Elapsed time is 17.875650 seconds.（节省时间，初始步长设为2.5）
    
13. >> tic,[mx,minf]=randwalk(z,[unifrnd(-1,7),unifrnd(-2,2)],2.5,0.000001,100,3),toc
    
14. mx =
    
15.    2.898269005888949  -0.857313207763948
    
16. minf =
    
17.   -0.023354083261446
    
18. Elapsed time is 11.866479 seconds.（节省时间，每步有3个位置供挑选）

复制代码

每步的随机挑选位置越多，步长越适中，全局寻优的能力就越强，当然，计算时间也越长，这些都在使用者可控范围内。大家有兴趣可以在1stOpt里试试，如果按算法默认设置，不自己选择算法，参数，看看寻优时间和结果。

[ 本帖最后由 rocwoods 于 2007-8-8 15:15 编辑 ]

shamohu

附件是一篇全局优化算法的文章：基于育种系统管理全局优化算法及其性能

有兴趣的试试，看看效果怎样！

qingruan

用以前写的PSO算法的程序跑了一下，基本上20次迭代内都能找到全局最优。程序就不贴了，用labview写的，以前在labview版贴过。把结果贴在这里，

顶上的黑点（变成jpg可能看不太清）就是找到的最优解。

rongrong-whf

现在用遗传和免疫解决函数优化问题的比较多，不知道楼主是否可以解释一下这种算法和您所说的算法之间的优劣？！

rocwoods

免疫算法我没有具体研究过。遗传算法我接触过，但是研究的不是很深入，就我个人的理解，遗传算法是一个非常巧妙的算法。这种巧妙体现在编码以及交叉算子的选择上。遗传算法和前面帖子里提到的随机行走法的一个共同特点是随机性，即程序在相同的初始条件下，每次运行完后的结果不一定是固定的，之所以每次结果不见得一样，是因为中间迭代过程中一些随机因素造成的。但是可以肯定的是，这两种算法都是向着全局最优的方向运行。
遗传算法是启发式算法，如果选择比较复杂的编码以及选择、交叉、变异算子的话，编程实现起来会有一定的难度。但是其适用范围非常广泛，一些离散的优化问题，利用传统优化算法很难求解的问题，利用遗传算法往往能收到不错的效果。
单纯说某一个算法优劣我觉得没多大意义，算法的好坏是和问题紧紧联系在一起的。每种算法都有其适合处理的问题。再强调一下，遗传算法的编码以及交叉算子的设计对算法的全局寻优能力以及收敛速度有很大影响。我们应该针对问题的特点设计编码方案以及选择、交叉、变异算子。

ljelly

我认为版主说的很有道理，遗传算法的确是一个不错的全局优化算法
无论是哪种算法，关键就是能够收敛，不会陷入局部极小，而达不到全局优化的目的
遗传算法我认为最重要的不是编码和算子的问题，最关键的应该是适应度函数的选取
决定了最终搜索的好坏，二进制编码要比十进制编码效果好，但麻烦，十进制则方便

rongrong-whf

感谢版主的解释，“我们应该针对问题的特点设计编码方案以及选择、交叉、变异算子”，这一点很重要，我深有体会。关键就是寻找适合自己问题的优化算法，但如何针对自己的问题选择算法呢？我接触这方面的比较少，还望版主指点！:D
    目前我在做简单的免疫算法，也是用于函数优化，只是函数较为复杂，刚好有现成的FORTRAN程序，所以想在MATLAB中调用一下，但遇到一些问题不知版主对于MATLAB调用FORTRAN程序了解否，可否指点？:P

wlk1978

本人也想在matlab中调用fortran进行优化,但有个问题没明白: matlab中用ga优化要写出函数形式的,现在我的函数要在fortran中经过一系列计算才能得到,而且其中要用到要优化的参数值,也就是要进行相互调用,不知道怎么实现?

wanglei5201118

我也遇到13楼的问题，我想这种问题还是很多见的啊，希望高手能给些建议！

mathcd

min=(@sin(@sqrt((x-50)^2+(y-50)^2)+@exp(1)))/(@sqrt((x-50)^2+(y-50)^2)+@exp(1))+1;
@bnd(0,x,100);
@bnd(0,y,100);
                       Variable           Value        Reduced Cost
                              X        51.74718            0.000000
                              Y        49.68623            0.000000


max=(@sin(@sqrt((x-50)^2+(y-50)^2)+@exp(1)))/(@sqrt((x-50)^2+(y-50)^2)+@exp(1))+1;
@bnd(0,x,100);
@bnd(0,y,100);
                       Variable           Value        Reduced Cost
                              X        50.00000            0.000000
                              Y        50.00000            0.000000

mathcd

这里提供三幅图，使用mathcd14绘出。

mathcd

版主提供的1stOpt帮助文档里给的针状函数，我在1stOpt1.5中很容易算出来，跟LINGO11差不多.但该软件宣称的“据测试，目前尚无其它软件，如著名的Lingo/Lindo，能求出此例的正确答案”的隐函数：
Minimum = z,z = sin((z*x-0.5)^2 + x*2*y^2-z/10)*exp(-((x-0.5-exp(-y+z))^2 + y^2-z/5+3))；
我用1stOpt1.5花了8分多钟还没有求出来，只好终止其运算。我下载了1stOpt1.5以及1stOpt的帮助文档，发现该软件的能力有夸大之嫌。对于普通的全局优化它并不比LINGO11好些，实际上它的计算时间比LINGO11要多，结果也没有那么神奇。其次，下图也是1stOpt帮助文档中的，我用LINGO11算过，用1stOpt1.5算过，用mathcd2001pro（初值完全是任意的）算过。从速度来说，1stOpt没有显示什么优势。从结果来看，优势在有无之间。

[ 本帖最后由 mathcd 于 2009-1-21 11:11 编辑 ]

rocwoods

1stOpt,我不熟，不好贸然下一个结论说好坏。就像有人MATLAB用不好就骂MATLAB哪儿哪儿不行，弱一样。软件是死的，人是活的。
每个软件有其特长之处，1stOpt优势就在于优化以及拟合等方面。楼上朋友可以找MathCAD Maple版本版主shamohu切磋1stOpt的使用问题，他比较在行。另外振动论坛dingd兄也比较熟悉。本版还有几个人比较熟悉它的使用，他们比我有发言权。
1stOpt也在发展，新版本比以前强很多，毕竟是国产的，还是要支持的。
楼上朋友可以用任何其他软件试下下面的问题：
http://forum.simwe.com/thread-806242-1-3.html

tangbaby

（我明明是回复的，怎么是新帖啊，请管理员移动到“全局最优化的讨论”）
LINGO求解，三行代码搞定了！求解耗时0.1MS
不知道我求的对不对我贴个截图（默认x>=0)

Ryan_P2011

楼主你好！
我在matlab中循环用fsolve解下列超越方程组，
以下是超越方程组，其中x1min和x2min是要求解的变量，Rou为已知量，其余的量如Proz4、a12、Matrix1等都是解一次方程就要变一次的量，Pro2从0~360变化，求出每个Pro值对应的x1min和x2min：

eq1=((cos(x2min)+1/2*sin(Proz4+x1min)^2*(1-cos(x2min)))*cos(x1min)+(-1/2*sin(Proz4+x1min)*cos(Proz4+x1min)*(1-cos(x2min))-1/2*2^(1/2)*sin(x2min))*sin(x1min))*Matrix1+(-(cos(x2min)+1/2*sin(Proz4+x1min)^2*(1-cos(x2min)))*sin(x1min)+(-1/2*sin(Proz4+x1min)*cos(Proz4+x1min)*(1-cos(x2min))-1/2*2^(1/2)*sin(x2min))*cos(x1min))*Matrix2+(1/2*sin(Proz4+x1min)*(1-cos(x2min))-1/2*2^(1/2)*cos(Proz4+x1min)*sin(x2min))*Matrix3-sin(Rou)*sin(Pro2)*cos(a12)-sin(Rou)*cos(Pro2)*cos(a22)+cos(Rou)*cos(a32);

eq2=((-1/2*sin(Proz4+x1min)*cos(Proz4+x1min)*(1-cos(x2min))+1/2*2^(1/2)*sin(x2min))*cos(x1min)+(cos(x2min)+1/2*cos(Proz4+x1min)^2*(1-cos(x2min)))*sin(x1min))*Matrix1+(-(-1/2*sin(Proz4+x1min)*cos(Proz4+x1min)*(1-cos(x2min))+1/2*2^(1/2)*sin(x2min))*sin(x1min)+(cos(x2min)+1/2*cos(Proz4+x1min)^2*(1-cos(x2min)))*cos(x1min))*Matrix2+(-1/2*cos(Proz4+x1min)*(1-cos(x2min))-1/2*2^(1/2)*sin(Proz4+x1min)*sin(x2min))*Matrix3-sin(Rou)*sin(Pro2)*cos(b12)-sin(Rou)*cos(Pro2)*cos(b22)+cos(Rou)*cos(b32);
前有用fsolve求解过，但是得出来的解反带回原方程组只有前几个点是正确的，不知道是否与初值的设置有关？如果是，要怎么解决？matlab还有其他解超越方程组的函数吗？
请楼主指教！先谢谢啦