非线性约束优化等式约束问题(fmincon)

laihai718

求一个最优化函数

max S(P)=-Σpi*ln(pi)
s.t. Σi*pi=E1      （ i=1：8）
     Σpi*(i-E1)^j=Ej   （ j=2：m）
     Σpi=1
pi>0
Ej是常量 代表第j阶矩的期望值
由样本得到，样本是4，2，6，5，7，6，7，2
E1=4.875，E2=27.375，E3=165.375，E4=1038.4


我是这样做的：
首先建立个m文件存放目标函数命名为myfun2
function f=myfun2(x)
f=x(1).*log(x(1))+x(2).*log(x(2))+x(3).*log(x(3))+x(4).*log(x(4))+x(5).*log(x(5))+x(6).*log(x(6))+x(7).*log(x(7))+x(8).*log(x(8));
保存后，在工作区一次输入：
>> A=[1 2 3 4 5 6 7 8;
    (1-4.875)^2 (2-4.875)^2 (3-4.875)^2 (4-4.875)^2 (5-4.875)^2 (6-4.875)^2 (7-4.875)^2 (8-4.875)^2;
    (1-4.875)^3 (2-4.875)^3 (3-4.875)^3 (4-4.875)^3 (5-4.875)^3 (6-4.875)^3 (7-4.875)^3 (8-4.875)^3;
    (1-4.875)^4 (2-4.875)^4 (3-4.875)^4 (4-4.875)^4 (5-4.875)^4 (6-4.875)^4 (7-4.875)^4 (8-4.875)^4;
    1 1 1 1 1 1 1 1];
>> b=[4.875 27.375 165.375 1038.4 1]';
>> p0=[0.15 0.1 0.15 0.1 0.1 0.15 0.15 0.1]';
>> [x,fval]=fmincon(@myfun2,p0,A,b)


我算的的结果都是0.125，答案是错的，不知道错在哪里了，请帮帮忙了，谢谢。

messenger

你的约束是等式约束，但[x,fval]=fmincon(@myfun2,p0,A,b) 是不等式约束，等式约束应为[x,fval]=fmincon(@myfun2,p0,[],[],A,b)

laihai718

感谢messenger 版主的热心解答，我运行了出现下面的警告，而且，程序一直在运行，一直都是busy

Warning: Trust-region-reflective method does not currently solve this type of problem,
using active-set (line search) instead.

最近也看了一些优化方面的资料，感觉就那几个命令，可怎么用就是不对，自己闷着头做了两三天了，好不容易算出一个结果吧，发现是错的，恼火

laihai718

我终止后，红字提示：
> In fmincon at 422
??? Operation terminated by user during ==> nlconst at 665

In ==> fmincon at 696
    [X,FVAL,LAMBDA,EXITFLAG,OUTPUT,GRAD,HESSIAN]=...

messenger

这个警告倒没什么，是说默认的求解方法不适合你的求解问题，求解器换另外一种方法求解，是求解器内部的事。

你的问题等式约束太多了，求不出结果来很正常，等式越多，能满足等式的结果就越少。

laihai718

我也觉得是等式太多了，我看的例子等式约束优化的比较少，即使有，也就一个等式约束的例子，晕，那这个题难不成matlab做不出来了？呵呵，还是感谢msg版主哈

wanglu

楼主这个问题我想用Forcal试一下，但我看不懂matlab代码，楼主能否用一般的数学语言描述一下？

wanglei5201118

看的不是很明白，建议用数学语言和自然语言重新描述一次！！这样大家才好帮你啊

messenger

其实，lz 最前面的描述就是数学语言描述，只是以比较简捷的形式描述的，下面以展开的形式再描述一下。



目标函数：min f(x)=x1*ln(x1)+x2*ln(x2)+x3*ln(x3)+x4*ln(x4)+x5*ln(x5)+x6*ln(x6)+x7*ln(x7)+x8*ln(x8)

约束函数：

x1+2*x2+3*x3+4*x4+5*x5+6*x6+7*x7+8*x8=E1;

x1*(1-E1)^2+x2*(2-E1)^2+x3*(3-E1)^2+x4*(4-E1)^2+x5*(5-E1)^2+x6*(6-E1)^2+x7*(7-E1)^2+x8*(8-E1)^2=E2;

x1*(1-E1)^3+x2*(2-E1)^3+x3*(3-E1)^3+x4*(4-E1)^3+x5*(5-E1)^3+x6*(6-E1)^3+x7*(7-E1)^3+x8*(8-E1)^3=E3;

x1*(1-E1)^4+x2*(2-E1)^4+x3*(3-E1)^4+x4*(4-E1)^4+x5*(5-E1)^4+x6*(6-E1)^4+x7*(7-E1)^4+x8*(8-E1)^4=E4;

x1+x2+x3+x4+x5+x6+x7+x8=1;

并且，x1>0;x2>0;x3>0;x4>0;x5>0;x6>0;x7>0;x8>0

其中，E1=4.875; E2=27.375; E3=165.375; E4=1038.4

laihai718

辛苦msg版主了，曾经用1stpot做过，但每次结果都不一样，很不稳定。

感谢大家一片热心。

wanglu

谢谢messenger 版主的问题描述。
我的Forcal解法是这样的：将x1,x2,x3作为优化参数，则由5个等式解线性方程组得到x4,x5,x6,x7,x8，用于计算目标函数。在约束函数中，也解线性方程组得到x4,x5,x6,x7,x8，验证其是否均大于0。优化采用自己的优化函数ROptMin，输出优化参数及5个等式的值。

代码：

1. 
   
2. x4_8(x1,x2,x3,x4,x5,x6,x7,x8:i:Ax55,Ax5)={  //求解线性方程组，x1,x2,x3为输入参数，x4,x5,x6,x7,x8为输出参数
   
3.   //E1=4.875, E2=27.375, E3=165.375, E4=1038.4,
   
4.   Ax55.XSLSF::setra[0, 4, 5, 6, 7, 8,
   
5.     (4-4.875)^2, (5-4.875)^2, (6-4.875)^2, (7-4.875)^2, (8-4.875)^2,
   
6.     (4-4.875)^3, (5-4.875)^3, (6-4.875)^3, (7-4.875)^3, (8-4.875)^3,
   
7.     (4-4.875)^4, (5-4.875)^4, (6-4.875)^4, (7-4.875)^4, (8-4.875)^4,
   
8.     1, 1, 1, 1, 1],
   
9.   Ax5.XSLSF::setra{0, 4.875-[x1+2*x2+3*x3], 27.375-[x1*(1-4.875)^2+x2*(2-4.875)^2+x3*(3-4.875)^2], 165.375-[x1*(1-4.875)^3+x2*(2-4.875)^3+x3*(3-4.875)^3], 1038.4-[x1*(1-4.875)^4+x2*(2-4.875)^4+x3*(3-4.875)^4], 1-[x1+x2+x3]},
   
10.   i=XSLSF::egaus(Ax55,Ax5),
    
11.   Ax5.XSLSF::getra{0,&x4,&x5,&x6,&x7,&x8},i
    
12. };
    
13. f(x1,x2,x3:x4,x5,x6,x7,x8:Ax55,Ax5)=    //目标函数定义
    
14. {
    
15.   x4_8(x1,x2,x3,&x4,&x5,&x6,&x7,&x8),
    
16.   x1*ln(x1)+x2*ln(x2)+x3*ln(x3)+x4*ln(x4)+x5*ln(x5)+x6*ln(x6)+x7*ln(x7)+x8*ln(x8)
    
17. };
    
18. s(x1,x2,x3:x4,x5,x6,x7,x8,i:Ax55,Ax5)=   //约束条件定义
    
19. {
    
20.   i=x4_8(x1,x2,x3,&x4,&x5,&x6,&x7,&x8),
    
21.   if{x4>0 & x5>0 & x6>0 & x7>0 & x8>0 & i, return(1)},0
    
22. };
    
23. main(:k,d,x1,x2,x3,x4,x5,x6,x7,x8, E1,E2,E3,E4:Ax55,Ax5)=
    
24. {
    
25.     E1=4.875, E2=27.375, E3=165.375, E4=1038.4,
    
26.     Ax55=new[rtoi(real_s),rtoi(5),rtoi(5)],
    
27.     Ax5=new[rtoi(real_s),rtoi(5)],
    
28.     k=50,x1=0.5,x2=0.3,x3=0.06,  //初值
    
29.     d=fcopt::ROptMin[HFor("f"),HFor("s"),10000,&k,1e-9,0,1: &x1,1e-300,1, &x2,1e-300,1 &x3,1e-300,1],
    
30.     x4_8(x1,x2,x3,&x4,&x5,&x6,&x7,&x8),
    
31.     printff{"\r\nx1={1,r}, x2={2,r}, x3={3,r}, x4={4,r}, x5={5,r}, x6={6,r}, x7={7,r}, x8={8,r}, f={9,r}, k={10,i}\r\nf0={11,r}, f1={12,r}, f2={13,r}, f3={14,r}, f4={15,r}, f5={16,r}\r\n",
    
32.           x1,x2,x3,x4,x5,x6,x7,x8,d,k,
    
33.           x1*ln(x1)+x2*ln(x2)+x3*ln(x3)+x4*ln(x4)+x5*ln(x5)+x6*ln(x6)+x7*ln(x7)+x8*ln(x8),
    
34.           [x1+2*x2+3*x3+4*x4+5*x5+6*x6+7*x7+8*x8-E1],
    
35.           [x1*(1-E1)^2+x2*(2-E1)^2+x3*(3-E1)^2+x4*(4-E1)^2+x5*(5-E1)^2+x6*(6-E1)^2+x7*(7-E1)^2+x8*(8-E1)^2-E2],
    
36.           [x1*(1-E1)^3+x2*(2-E1)^3+x3*(3-E1)^3+x4*(4-E1)^3+x5*(5-E1)^3+x6*(6-E1)^3+x7*(7-E1)^3+x8*(8-E1)^3-E3],
    
37.           [x1*(1-E1)^4+x2*(2-E1)^4+x3*(3-E1)^4+x4*(4-E1)^4+x5*(5-E1)^4+x6*(6-E1)^4+x7*(7-E1)^4+x8*(8-E1)^4-E4],
    
38.           [x1+x2+x3+x4+x5+x6+x7+x8-1]
    
39.           },
    
40.     delete[Ax55],delete[Ax5]
    
41. };
    

复制代码

问题：经多次尝试验证，无法获得一组符合约束条件的x1,x2,x3,x4,x5,x6,x7,x8，故求解无法进行。
楼主若能提供一组符合约束条件的x1,x2,x3,x4,x5,x6,x7,x8，估计优化可以进行下去。

laihai718

>> p0=[0.15 0.1 0.15 0.1 0.1 0.15 0.15 0.1]'

这组就符合条件哦

wanglu

x1=0.14999999999999999, x2=0.10000000000000001, x3=0.14999999999999999, x4=6.2367289225260389, x5=-30.928361002604152, x6=66.222521972656239, x7=-63.881876627604164, x8=22.950986735026042

这是取x1=0.14999999999999999, x2=0.10000000000000001, x3=0.14999999999999999所求的x4,x5,x6,x7,x8，不符合约束x1>0;x2>0;x3>0;x4>0;x5>0;x6>0;x7>0;x8>0

如果p0=[0.15 0.1 0.15 0.1 0.1 0.15 0.15 0.1]'，则其余4个等式可能不成立（没有验证）。

我说的初值要满足所有约束，似乎很难得到？

laihai718

感谢您哦，我有到是有一个，是准确的值，如果贴出来了，就没意义了，您所要的初值，我还真不好给您弄来。

messenger

总感觉你的哪个约束条件写错了，我算了几次总是提示说约束条件冲突。

你确定你的准确值可以满足这5个约束条件？

要不然你把准确值贴出来吧，这样便於大家调试程序，反正你留着准确值也没什么用，如果准确值有用，你也就不用再费事算了。

laihai718

p1=0.09135128
p2=0.09211282
p3=0.10783846
p4=0.13826154
p5=0.17520513
p6=0.18946152
p7=0.14442821
p8=0.06134103

这个就是准确值，辛苦大家了，在此感谢

wanglu

楼主这组解带入5个等式，只要最后一个等式成立，难道我算错了？
f1=-0.20834106000000041, f2=-23.263215382656249, f3=-170.15476392861328, f4=-1000.9975542669556, f5=-9.9999998282029878e-009

laihai718

感谢楼上的这么热心，我的约束条件是
Σpi=1
pi>0

您解的好像不满足这两个条件吧

messenger

我怎么觉得你的问题本身就没解呢。17#已经说了，你的精确值带入等式，也不满足等式呀。

laihai718

你们说的这些个问题，我还是验证了才知道的，之前没想着答案会有问题，这个是我从一本书上看到的，估计是印刷错误，我还得好好研究研究，真是麻烦你们了，谢谢各位。