1stOpt一个例子的说明求解释，另结果是否最优？

wanglu

首先，这个例子似乎最早从网上搜来，什么地方不记得了，故放到了Forcal优化库FcOpt的中【3.1.3 选自网上的难题例子】：http://www.forcal.net/sysm/forcal9/fchtm/fcopt.htm

在FcOpt本次更新时，用新增函数Opt算出的结果比以前更好，故更新了这个例子。
一直想再找找这个例子，看看以前的结果，却偶然在1stOpt中发现了，不管是哪里的例子，先表示感谢！

该例子在1stOpt使用手册2009，【2.5.4 关键字“PassParameter”的使用】
例 非线性约束拟合
拟合数据
x -0.08,-0.065,-0.05,-0.03,-0.015,0.015,0.03,0.05,0.065,0.08,0
y 20.26008,19.72613,19.501619,18.72662,18.58769,18.592199,18.88372,19.5453,19.88743,20.9914,18.12336
拟合公式
y = a * (x − d )^4 + b * (x − d)^2 + c

其中：x、y：自变量和因变量；a、b、c、d：待求参数。

约束条件：
1）必须大于0；
2）计算y 与实际y 的最大差值的绝对值不大于0.25，即：Max(abs( yi - yi' )) <=0.25, i =1..11

有一处说明，感到不解：1）必须大于0；但所给结果有些却小于0，为什么？

另外，感到不加约束时1stOpt的结果并非最优，1stOpt的结果：
目标函数值(最小): 0.309952381552963
a: -6877.78414257241
b: 385.898376776761
c: 18.4260680489247
d: -0.00393418836172759

使用与1stOpt相同的目标函数，Forcal的最优结果：

1. !using["fcopt","math"];
   
2. init(::Array,max)=
   
3. {
   
4.     max=11,
   
5.     Array=arrayinitns{max,2 :
   
6. "-0.08 20.26008
   
7. -0.065 19.72613
   
8. -0.05 19.501619
   
9. -0.03 18.72662
   
10. -0.015 18.58769
    
11. 0.015 18.592199
    
12. 0.03 18.88372
    
13. 0.05 19.5453
    
14. 0.065 19.88743
    
15. 0.08 20.9914
    
16. 0 18.12336"
    
17.     }.free()
    
18. };
    
19. f(a, b, c, d :i,s,x,y:Array,max)=
    
20. {
    
21.     s=0,i=0,(i<max).while{
    
22.         x=Array[i,0], y=Array[i,1],
    
23.         s=s+[ a*(x-d)^4+b*(x-d)^2+c-y]^2,
    
24.         i++
    
25.     },
    
26.     s
    
27. };
    
28. Opt[HFor("f"), optdeep,50];

复制代码

结果：
420.1374351331109         -170.7494237822797        35.81332551003671         -0.4495641127890136       0.2619659051106249

如果确实如此，希望1stOpt能更新此例子。当然，一个例子说明及结果的偶然疏漏，不能丝毫减弱1stOpt在优化领域的领先地位，任何东西都是瑕不掩瑜的。或者这个例子我还没有看懂，望解释。

wanglu

这个网页中有相同的例子：http://cos.name/cn/topic/7605
但似乎还不是我最初看到的那个。

shamohu

看了一下原手册，代码中有一句“Parameter a=[,0],b,c,d;”，也即参数a必须小于0，估计是手册文本笔误，导致与代码不一致，代码如果改为“Parameter a,b,c,d;”，检验代码如下：

1. 
   
2. Parameter a,b,c,d;
   
3. Function y = a*(x-d)^4+b*(x-d)^2+c;
   
4. Data;
   
5. x=-0.08,-0.065,-0.05,-0.03,-0.015,0.015,0.03,0.05,0.065,0.08,0;
   
6. y=20.26008,19.72613,19.501619,18.72662,18.58769,18.592199,18.88372,19.5453,19.88743,20.9914,18.12336;
   

复制代码

100%概率得到结果：

均方差(RMSE): 0.154321360319194
残差平方和(SSE): 0.261965904758431
相关系数(R): 0.981842657141946
相关系数之平方(R^2): 0.964015003383557

参数        最佳估算
----------        -------------
a        420.080615896829
b        -170.748499130253
c        35.8154870740463
d        -0.449593738487319

运行了一下1楼的Forcal代码，10次结果：
1: 420.5617281887409         -170.7523819414472        35.79641084920254         -0.4493379602386528       0.261965925167395
2: -6877.911075436947        385.8994487490645         18.42606766813            -3.934273670173601e-003   0.3099523816639323
3: -6877.911075436947        385.8994487490645         18.42606766813            -3.934273670173601e-003   0.3099523816639323
4: -6877.514585672865        385.8960688679966         18.42607125320205         -3.934125343761715e-003   0.3099523816443295
5: 420.110108154518          -170.7492626641425        35.81442052315673         -0.4495786747918449       0.2619659049557511
6: 420.5617281887409         -170.7523819414472        35.79641084920254         -0.4493379602386528       0.261965925167395
7: -6877.45424461292         385.8959942050641         18.42607459032513         -3.934109677197721e-003   0.3099523818516957
8: -6878.624930337377        385.9032095700601         18.42606421519121         -3.934076759786656e-003   0.3099523819268166
9: 420.5617281887409         -170.7523819414472        35.79641084920254         -0.4493379602386528       0.261965925167395
10: -6878.212494911696        385.9010986236473         18.42606836119495         -3.934202499195048e-003   0.309952381648598

也即大概50%的成功概率，已经不错了！

shamohu

另外，将公式改为：y = a*(x-d)^4+b*(x-d)^p2+c，还是用一楼的Forcal代码，结果大概是：
[a,b,c,d,p2] = 48591.8741934397          -6.374606769898301e-017   18.75746239720849         -3.234708294999146e-003   7.524779375888238e+017

从理论上，为使“*(x-d)^p2”项能进行计算，因为p2是非整实数，“x-d”必须为正才行，但Forcal得出的d为-3.234708294999146e-003，有几个点使得“x-d”为负，而计算还能照常进行，不知什么原理？是自动按复数还是其它？

wanglu

另外，将公式改为：y = a*(x-d)^4+b*(x-d)^p2+c，还是用一楼的Forcal代码，结果大概是：
[a,b,c,d,p2] = 48591.8741934397          -6.374606769898301e-017   18.75746239720849         -3.234708294999146e-003 ...
shamohu 发表于 2011-2-15 15:37

指数计算我用的是C++的pow函数，是不是p2=7.524779375888238e+017，指数太大了，pow就返回一个默认值，验证了一下：
(-0.08+3.234708294999146e-003)^7.524779375888238e+017 ＝ 0

增加了p2后，拟合难度增加了，故修改代码如下：

1. !using["fcopt","math"];
   
2. init(::Array,max)=
   
3. {
   
4.     max=11,
   
5.     Array=arrayinitns{max,2 :
   
6. "-0.08 20.26008
   
7. -0.065 19.72613
   
8. -0.05 19.501619
   
9. -0.03 18.72662
   
10. -0.015 18.58769
    
11. 0.015 18.592199
    
12. 0.03 18.88372
    
13. 0.05 19.5453
    
14. 0.065 19.88743
    
15. 0.08 20.9914
    
16. 0 18.12336"
    
17.     }.free()
    
18. };
    
19. f(a, b, c, d, p2 :i,s,x,y:Array,max)=
    
20. {
    
21.     s=0,i=0,(i<max).while{
    
22.         x=Array[i,0], y=Array[i,1],
    
23.         s=s+[ a*(x-d)^4+b*(x-d)^p2+c-y]^2,
    
24.         i++
    
25.     },
    
26.     s
    
27. };
    
28. Opt[HFor("f"), optmax,300, optwaysimdeep, optwayconfra];

复制代码

运行几次后可得解：
421.4346987551445         -168.2006461102398        35.76126549040257         -0.4467645916460592       1.981167282514616         0.2619657217171444

若使用以下代码，缩小搜索范围，可以一半以上的概率获得结果：

1. !using["fcopt","math"];
   
2. init(::Array,max)=
   
3. {
   
4.     max=11,
   
5.     Array=arrayinitns{max,2 :
   
6. "-0.08 20.26008
   
7. -0.065 19.72613
   
8. -0.05 19.501619
   
9. -0.03 18.72662
   
10. -0.015 18.58769
    
11. 0.015 18.592199
    
12. 0.03 18.88372
    
13. 0.05 19.5453
    
14. 0.065 19.88743
    
15. 0.08 20.9914
    
16. 0 18.12336"
    
17.     }.free()
    
18. };
    
19. f(a, b, c, d, p2 :i,s,x,y:Array,max)=
    
20. {
    
21.     s=0,i=0,(i<max).while{
    
22.         x=Array[i,0], y=Array[i,1],
    
23.         s=s+[ a*(x-d)^4+b*(x-d)^p2+c-y]^2,
    
24.         i++
    
25.     },
    
26.     s
    
27. };
    
28. Opt[HFor("f"), optwaysimdeep, optwayconfra, optrange,-1e10,1e10,-1e10,1e10,-1e10,1e10,-1e10,1e10,-1e10,1e10];

复制代码

结果：
421.4345728315115         -168.2007599061191        35.7612703620313          -0.4467647520367266       1.981168132343698         0.2619657217171494

不知1stOpt最优结果如何？

shamohu

1stOpt能以100%概率得到最优解，用时也很短。

Forcal似乎很深奥，各种控制参数不好懂也不好掌握，能否简化直观些？

5#的两段代码计算都比较耗时，成功率也不高。

“(-0.08+3.234708294999146e-003)^7.524779375888238e+017 ”在C++中结果为0，估计是把“7.524779375888238e+017 ”自动视为整数了，用其它编译器会报错。

wanglu

1stOpt能以100%概率得到最优解，用时也很短。

Forcal似乎很深奥，各种控制参数不好懂也不好掌握，能否简化直观些？

5#的两段代码计算都比较耗时，成功率也不高。

“(-0.08+3.234708294999146e-003)^7.52477 ...
shamohu 发表于 2011-2-15 16:54

Forcal不是专门做优化的，只是喜欢这个就写了这个库。如果像1stOpt那样，给Opt函数一个界面，参数选择就方便了。实际上大多数参数不需要修改，我经常改的就是上面出现的几个参数：optmax、optwaysimdeep、optwayconfra、optrange。

也想设计成1stOpt的模式，给出公式、数据就可求解，但想想太麻烦，又失去了灵活性。
感觉用户自己编写目标函数虽然稍稍麻烦了些，但非常灵活，也不需要添加那么多的关键字（关键字多了也会增加复杂性的），什么含积分的优化、微分方程优化、复数方程优化等等都可以适应。

仍然上例，如果公式改成：y=a*(x-d)^p1+b*(x-d)^p2+c ，不知1stOpt结果如何？
Forcal这个结果是否最优：
3087.197326547206         -2736.083082665887        20.26150232509316         -8.000000000000457e-002   1.666935681811038         1.603644423125645         0.2283370225460334

wanglu

另外，shamohu 版主是否用过Fortran，我有个问题可否帮忙看一下（或者其他人或其他人的朋友帮忙看一下）：http://forum.simwe.com/thread-969615-1-1.html

shamohu

均方差(RMSE): 0.144073533164399
残差平方和(SSE): 0.228329012543204
相关系数(R): 0.984192845971347
相关系数之平方(R^2): 0.968635558061179

参数        最佳估算
----------        -------------
a        3751.74907707895
b        -3400.70265100457
c        20.2611557025148
d        -0.080000245145933
p1        1.66107056588994
p2        1.60957322264913

将公式改为y=a*(x-d)^p1+b*(x-d)^p2+c，运行5#的两段代码，均给出错误信息：
Forcal表达式运行错误！
实数表达式名称：  所在模块：-53
实数函数名：fcopt::Opt  错误代码：4

wanglu

因增加了参数p1，需做如下修改：

1. f(a, b, c, d, p1, p2 :i,s,x,y:Array,max)=
   
2. {
   
3.     s=0,i=0,(i<max).while{
   
4.         x=Array[i,0], y=Array[i,1],
   
5.         s=s+[ a*(x-d)^p1+b*(x-d)^p2+c-y]^2,
   
6.         i++
   
7.     },
   
8.     s
   
9. };
   
10. Opt[HFor("f"), optwaysimdeep, optwayconfra, optrange,-1e10,1e10,-1e10,1e10,-1e10,1e10,-1e10,1e10,-1e10,1e10,-1e10,1e10];

复制代码

不过得到1stOpt的结果不容易。