失之毫厘谬以千里的优化测试题！

shamohu

下面这道优化测试题看似简单，但想得到最优解却也不易，较之国际上经典的测试题诸如Griewangk's function，Schaffer Function等难多了，有兴趣的可先试一下除1stOpt以外的任何其它软件或自己编程，看能得到什么样的结果。

min ＝ (x1^x2+x2^x1-5*x1*x2*x3-85)^2+(x1^3*x2^x3*x3^x2-60)^2+(x1^x3+x3^x1-x2-0.55)^2

rocwoods

我算出来的结果是

1. x1= 5.176045131038809,
   
2. x2=2.301301559119753,
   
3. x3=0.566071813266631  minf=0.0690
   
4. x1=4.837351077980695
   
5. x2=2.438534929133349
   
6. x3=0.615526123391624 minf=0.0646
   
7. 应该有多组极小值，方法：先用蒙特卡罗法找到一个最小点，然后用该点作为初始值，用我以前发的randwalk函数或者MATLAB自带的fminunc来求解（fmincun的option设置下可能得到更好的结果）。说明，有时候不成功，busy时间长了，ctrl+c中断，多试两回。
   
8. 代码如下：
   
9. f=@(x)  (x(1).^x(2)+x(2).^x(1)-5*x(1).*x(2).*x(3)-85).^2+(x(1).^3.*x(2).^x(3).*x(3).^x(2)-60).^2+(x(1).^x(3)+x(3).^x(1)-x(2)-0.55).^2;f
   
10. x=unifrnd(0,10,100000,3);
    
11. F=zeros(100000,1);
    
12. for k=1:100000
    
13. F(k)=f(x(k,:));
    
14. end
    
15. [Fmin,ind]=min(F);
    
16. minx=x(ind,:);
    
17. [xmin,minf]=fminunc(f,minx)
    
18. 或者
    
19. f=@(x)  (x(1).^x(2)+x(2).^x(1)-5*x(1).*x(2).*x(3)-85).^2+(x(1).^3.*x(2).^x(3).*x(3).^x(2)-60).^2+(x(1).^x(3)+x(3).^x(1)-x(2)-0.55).^2;f
    
20. x=unifrnd(0,10,100000,3);
    
21. F=zeros(100000,1);
    
22. for k=1:100000
    
23. F(k)=f(x(k,:));
    
24. end
    
25. [Fmin,ind]=min(F);
    
26. minx=x(ind,:);
    
27. [mx,minf]=randwalk(f,minx,0.1,0.0001,100,3,3)
    
28. randwalk的代码如下：
    
29. function [mx,minf]=randwalk(f,x,lamda,epsilon,N,n,r)
    
30. %随机行走法求函数的极小值。输入f为所求函数的句柄，
    
31. %x为初始值。lamda为步长。epsilon为控制lamda的减小的阈值，即lamda减小到epsilon时迭代停止。
    
32. %N为为了产生较好点的迭代控制次数。
    
33. %n为单步循环行走次数，目的是尽可能走到全局最优点附近
    
34. %函数返回值mx为n次试验得到的最优解，minf为响应的最优值。
    
35. %r,函数的维数
    
36. F=zeros(n,1);
    
37. X=zeros(n,r);
    
38. f1=f(x);
    
39. while(lamda>=epsilon)
    
40.     k=1;
    
41.     while(k<=N)
    
42.         u=unifrnd(-1,1,n,r);
    
43.         for ii=1:n
    
44.         X(ii,:)=x+lamda*(u(ii,:)/norm(u(ii,:)));
    
45.         F(ii)=f(X(ii,:));
    
46.         end
    
47.         [f11,kk]=min(F);
    
48.         if f11<f1 %&& X(kk,1)>=13 && X(kk,1)<=100 && X(kk,2)>=0 && X(kk,2)<=100 && ((X(kk,1)-5)^2+(X(kk,2)-5)^2>=100)
    
49.             f1=f11;
    
50.             x=X(kk,:);
    
51.             k=1;
    
52.         else
    
53.             k=k+1;
    
54.         end
    
55.     end
    
56.     lamda=lamda/2;
    
57. end
    
58. mx=X(kk,:);
    
59. minf=f1;

复制代码

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

shamohu

先谢过rocwoods版主的解答，不过离正解还有相当距离。

因为Lingo不需初值，试了一下Lingo，不论是否用选项“Use Global Solver”，计算的结果和rocwoods的都差不多：

Global optimal solution found.
  Objective value:                             0.6463502E-01
  Extended solver steps:                             341
  Total solver iterations:                        175730


                       Variable           Value        Reduced Cost
                             X1        4.834151           0.1967419E-07
                             X2        2.439950           0.5887098E-07
                             X3       0.6160029          -0.4238520E-07

但都不是正解啊！大家再试试！看下标题，应该知道结果不会太一般。

rocwoods

的确很不一般，呵呵。从求出的最优解来看，这个优化问题似乎等价于求解一个三元非线性方程组，我把它改成三元非线性方程组后，求解了一番，得到的结果还是上面我给出的。后来想到用solve符号求解试试，算了一晚上，早晨睡觉醒来居然还在算。既没有提示无解也没有给出近似解。
看来得寻找有效的算法。

hnyzyhq230

这个优化问题的x1、x2、x3有上下限范围吗?

shamohu

没有范围限制。

smarten

我不懂优化，看看Mathematica, 算的也是0.064635。

1. f = (x1^x2 + x2^x1 - 5*x1*x2*x3 - 85)^2 + (x1^3*x2^x3*x3^x2 -
   
2.      60)^2 + (x1^x3 + x3^x1 - x2 - 0.55)^2

复制代码

1. NMinimize[{f, x1 > 0, x2 > 0, x3 > 0}, {x1, x2, x3},
   
2. Method -> "RandomSearch"]

复制代码

得到：

1. {0.064635, {x1 -> 4.83415, x2 -> 2.43995, x3 -> 0.616003}}

复制代码

或者：{0.064635, {x1 -> 4.83415, x2 -> 2.43995, x3 -> 0.616003}}
这个题目关键是加上x1 > 0, x2 > 0, x3 > 0的限制之后，Mathematica的那些方法都用不上，"LevenbergMarquardt"方法应该更好一些。
可能需要一些变通。

bainhome

Mathematica和MATLAB的结果很多时候是差不多的，LINGO有的时候优化解能够更好，但都不如1stopt的随机初值算法凶猛。
这么长时间，shamohu兄可以给答案了，谢谢！

shamohu

2.0以前的版本不敢担保，但2.0, 2.5版运行下面代码，大概80-90％的概率能得到最优解。

Algorithm = GLM[200];
MinFunction (x1^x2+x2^x1-5*x1*x2*x3-85)^2+(x1^3*x2^x3*x3^x2-60)^2+(x1^x3+x3^x1-x2-0.55)^2;

最优结果：
目标函数值(最小): 0
x1: 19390.8841948436
x2: 0.45
x3: 2.33057529481628E-25

如果将x3变为2.33057529481628E-24，从绝对值来看变化是很微小的了，但目标函数变为11903.45934，这也是题目所言“失之毫厘，谬以千里”之由。

下面是Lingo代码，注意初值已经给的非常接近1stOpt得到的最优值了，但最终结果还是回到了局部最优（11.0试用版）。

Init:
x1= 19390;
x2= 0.45;
x3= 2.33057529481628E-25;
EndInit
Min =(x1^x2+x2^x1-5*x1*x2*x3-85)^2+(x1^3*x2^x3*x3^x2-60)^2+(x1^x3+x3^x1-x2-0.55)^2;
Global optimal solution found.
  Objective value:                             0.6463502E-01
  Objective bound:                             0.6463500E-01
  Infeasibilities:                              0.000000
  Extended solver steps:                             571
  Total solver iterations:                        210589

                       Variable           Value        Reduced Cost
                             X1        4.834151           0.1778444E-07
                             X2        2.439950           0.2640767E-07
                             X3       0.6160029           0.4128631E-07

这道题仅是一例，并不能由此下结论孰优孰劣，每个软件都有优点和缺点。

[ 本帖最后由 shamohu 于 2008-8-7 10:07 编辑 ]

mathcd

版主的题目是：min(x1^x2+x2^x1-5*x1*x2*x3-85)^2+(x1^3*x2^x3*x3^x2-60)^2+(x1^x3+x3^x1-x2-0.55)^2;
我用lingo11计算：
代码：
Init:
x1= 19390;
x2= 0.45;
x3= 2.33057529481628E-25;
EndInit
Min =(x1^x2+x2^x1-5*x1*x2*x3-85)^2+(x1^3*x2^x3*x3^x2-60)^2+(x1^x3+x3^x1-x2-0.55)^2;
结果如下：
Global optimal solution found.
  Objective value:                             0.6463502E-01
  Objective bound:                             0.6463500E-01
  Infeasibilities:                              0.000000
  Extended solver steps:                             571
  Total solver iterations:                        210589


                       Variable           Value        Reduced Cost
                             X1        4.834151           0.1778444E-07
                             X2        2.439950           0.2640767E-07
                             X3       0.6160029           0.4128631E-07

                            Row    Slack or Surplus      Dual Price
                              1       0.6463502E-01       -1.000000

我用1st0pt1.5计算：
代码：
Algorithm = GLM[200];
MinFunction (x1^x2+x2^x1-5*x1*x2*x3-85)^2+(x1^3*x2^x3*x3^x2-60)^2+(x1^x3+x3^x1-x2-0.55)^2;
结果如下：
====== 结果 ======

迭代数: 189
计算用时(时:分:秒:毫秒): 00:00:00:500
计算中止原因: 达到收敛判定标准
优化算法: 标准简面体爬山法 + 通用全局优化法
函数表达式: (x1^x2+x2^x1-5*x1*x2*x3-85)^2+(x1^3*x2^x3*x3^x2-60)^2+(x1^x3+x3^x1-x2-0.55)^2
目标函数值(最小): 0.064635017964009
x1: 4.83415105079868
x2: 2.439950006258
x3: 0.616002912636042

====== 计算结束 ======
我再用mathcd2001pro计算：


有兴趣的可以自己比较。

mathcd

版主的题目用lingo11算时，如果就x1，x2，x3均设为free，则lingo11很难给出结果。
代码：
代码：
Init:
x1= 19390;
x2= 0.45;
x3= 2.33057529481628E-25;
EndInit
Min =(x1^x2+x2^x1-5*x1*x2*x3-85)^2+(x1^3*x2^x3*x3^x2-60)^2+(x1^x3+x3^x1-x2-0.55)^2;
@free（x1）;@free(x2);@free(x3);
如图，lingo11算了差不多30分钟，还是不能给出结果。


而用mathcd2001pro算出的最小值是最小的，优于1st0pt1.5和lingo11。

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

shamohu

1：Lingo的最新版是11，1stOpt是2.5，比较时可能的话最好都用最新版
2：9楼给出的结果是：

最优结果：
目标函数值(最小): 0
x1: 19390.8841948436
x2: 0.45
x3: 2.33057529481628E-25

而且初值是随机给的。MathCAD如果不用1stOpt的结果作初值，不知结果如何？

mathcd

MathCAD2001pro如果不用1stOpt的结果作初值，结果很糟。

wanglu

看到这个题，忍不住试了一下，Forcal代码：

1. f(x1,x2,x3)=(x1^x2+x2^x1-5*x1*x2*x3-85)^2+(x1^3*x2^x3*x3^x2-60)^2+(x1^x3+x3^x1-x2-0.55)^2;
   
2. main(:f,x1,x2,x3)=
   
3. {
   
4.   f=fcopt::GOpt0[HFor("f") : &x1, &x2, &x3],
   
5.   printff{"\r\nx1={1,r}, x2={2,r}, x3={3,r}, f={4,r}\r\n", x1,x2,x3, f}
   
6. };

复制代码

最好的结果：
x1=7045.8526271895898, x2=0.50141735691227995, x3=3.3412071021255034e-020, f=2.644090963625534e-003

1stopt好强哦。

wanglu

新年第一帖，做个老题。
再次看到这个题，仍忍不住试了一下，Forcal代码：

1. f(x1,x2,x3)=(x1^x2+x2^x1-5*x1*x2*x3-85)^2+(x1^3*x2^x3*x3^x2-60)^2+(x1^x3+x3^x1-x2-0.55)^2;
   
2. fcopt::Opt[HFor("f")];

复制代码

结果：
19390.88419484345         0.4500000000000003        2.330575294816493e-025    1.109335647967048e-031

wujianjack2

   第一次回贴，心情有点激动，大家莫怪。
   shamohu版主大人的这个问题，确实很经典啊！我偶然在一篇名为"1stOpt and Global Optimization Platform-Comparsion and Case Study"的文献中看到了这个测试题。看了这么多大家的解答，我发现大家求出的结果，x1,x2,x3都是非负的，即使最初题目中没有限制x1,x2,x3的取值范围。那么我想，这个问题的结果是否就是全为正解呢？
   大家还有一个共同点就是，在使用LINGO编写代码时，采用了开门见山式的编写方式，求出的结果不是很理想，然而，基于我上面结果全为非负的假设，是否可以使用LINGO求出较好解呢？

如下代码回答了这个问题：（此方案非我原创，特此说明）
[C1] X1^X2 + X2^X1 - 5*X1*X2*X3 - 85 = 0;
X1 = @EXP(LX1); @FREE(LX1);
X2 = @EXP(LX2); @FREE(LX2);
X3 = @EXP(LX3); @FREE(LX3);
[LC2]  3*LX1 + X3*LX2 + X2*LX3 - @LOG(60) = 0;
[C3] X1^X3 + X3^X1 - X2 - 0.55 = 0 ;

无论是使用LINGO默认的Local Nonlinear Solver还是Global Solver均可迅速得到如下结果：（软件版本为LINGO 11，抱歉，我用的是Crack,此问题Demo版本无法用Global Solver）
  Feasible solution found.
  Infeasibilities:                             0.1268152E-05
  Total solver iterations:                            28

                                           Variable           Value
                                                 X1        19390.88
                                                 X2       0.4500000
                                                 X3        0.000000
                                                LX1        9.872558
                                                LX2      -0.7985077
                                                LX3       -56.71851

                                                Row    Slack or Surplus
                                                 C1        0.000000
                                                  2      -0.1268152E-05
                                                  3        0.000000
                                                  4        0.000000
                                                LC2        0.000000
                                                 C3        0.000000
（以上结果为默认的Local Solver求解结果）。
看来，适当进行变化并进行限制，以适应LINGO求解的习惯，LINGO也能较好地求出这个问题的解。
我是一名本科生，所知甚少，冒然回答，希望各位坛友给出意见与建议。

lipenggg

1.5版的计算结果：
====== 结果 ======

迭代数: 543
计算用时(时:分:秒:毫秒): 00:00:09:893
计算中止原因: 达到收敛判定标准
优化算法: 麦夸特法(Levenberg-Marquardt) + 通用全局优化法
函数表达式: (x1^x2+x2^x1-5*x1*x2*x3-85)^2+(x1^3*x2^x3*x3^x2-60)^2+(x1^x3+x3^x1-x2-0.55)^2
目标函数值(最小): 8.31195618737763E-26
x1: 19390.629582248
x2: 0.450000597661798
x3: 2.33095483355795E-25