变尺度法BFGS算法的C--源码

=73这个是我那时候交的作业（很粗了。。）。嗯。。。应该标注很清楚了。不过如果你用matlab编写的话，很多大段的代码，基本上一句话就代替了。;E*M3R3{!C(W/c//变尺度法BFGS算法的C++源码#includeiostream.h#includemath.h#defineA85V*k;O-[0f!s#defineB10#defineC75E-j+L4~6x!R#defineD0#defineSTART{5,4,5,5}//初始迭代点0i7C9K7G%f7Qconstintn=4;//变量个数6\9q1`0c3h&H:gconstintnc=3;//不等式约束个数(也可根据题目需要增加等式约束个数修改约束构造函数)3d9O$C#p'p#k'\9yconstdoubleepsilon1=1E-6;//一维搜索精度4[0L(?&N1V0y8\constdoubleepsilon2=1E-4;//变尺度法收敛精度constdoubleepsilon=1E-4;//惩罚函数法收敛精度,j/{t3^,a%p7q:Adoubler=1;//惩罚因子doublebr=1;//变尺度法强行退出循环因子(可去除)constdoublezr=10;//递增系数doubledet=0.0014;//外惩罚函数法约束函数余量//(此det值仅适合当前ABCD，修改请将其恢复为0.0001或其他建议值)//计算梯度voidcomput_grad(double(*pf)(double*x),intn,double*point,double*grad);'d,f^-q0b4}$t3L)t4X//无约束BFGS变尺度法（目标函数，维数，最小点）返回0值doubleBFGS(double(*pf)(double*x),intn,double*min_point);//构造混合外点惩罚函数返回惩罚函数值（构造点）9v/v*K(x#o4V~)[doublefunP(double*p);w5L2u'C7BS#e!M//返回目标函数值,并返回约束函数值存入ys（计算点，约束函数值存放处）doublefun(double*x,double*ys);)\.b#l9V;a%c#x/};q5B#w//一维搜索函数返回最优步长（目标函数，维数，初始点，方向）doubleline_search(double(*pf)(double*x),intn,double*start,double*direction);'k.B7Lhv-`;Y(W//配给与一维搜索函数返回下一点的的惩罚函数值（计算点，计算方向，步长）:K-?-v(L1o-ddoublefunPt(double*x,double*direction,doublet);2n)V:F9G%W#?/F-Z9pBvoidmain(){inti,j;doublec=0;//表示外惩罚函数循环次数doublef_min,f_temp,f_punish;0tT;E+M2l5n9yk7mdouble*ys=newdouble[nc];O/X%S3d6v!W.^doublex_x=0,f_ys=0;doublex_min[n]=START;doublex_temp[n]=START;2e*n$sL#Sf_temp=fun(x_temp,ys);&I8t3^'q;pfor(;;){c++;-z.[4~3^/M-hBFGS(funP,n,x_min);-P7w9{#p9i'r3Z#N/n+S5\f_min=fun(x_min,ys);for(j=0;jn;j++)l!Z+L!z%Sx_x=x_x+pow((x_min[j]-x_temp[j]),2);S9d$q$r*W$`#i$V'W,pfor(i=0;inc;i++)f_ys+=fabs(ys);x_x=sqrt(x_x);:r:P;n2z0P/M:{+G2s//参考值显示cout\n外惩罚函数循环迭代第c轮endl目标函数:f_min=f_minendl自变量:x1=x_min[0]x2=x_min[1]x3=x_min[2]x4=x_min[3]endl3B1P*e*n7v7[+D约束函数:g1=ys[0]+detg2=ys[1]+detg3=ys[2]+detendl5l0z2i6^(f%K&l+k%Q判断收敛因子:x_x=x_xf_f=fabs(f_min-f_temp)det_f=fabs((f_min-f_temp)/f_min)f_ys=f_ysendl;3]*d5Z-P;[!x2S0d#D4n//-------------------|9[:`7k#m+h&m3@,E//收敛判断if((fabs(x_x)=epsilon)|(f_ys=epsilon)){cout因为迭代xk与xk+1距离x_x或者约束函数值f_ys达到收敛精度epsilon而认为收敛endl;!o4e#N5w3Obreak;}if((fabs(f_min)=1)&&(fabs(f_min-f_temp)=epsilon))+J%B!W1yd:w4_$z.^{cout因为迭代xk与xk+1的函数值差距f_f达到收敛精度epsilon而认为收敛endl;break;}0b4@1I-C.Q,XD;M8F'Delseif((fabs(f_min)=1)&&(fabs((f_min-f_temp)/f_min)=epsilon)){cout\n因为迭代xk与xk+1的函数值差距det_f达到收敛精度epsilon而认为收敛endl;:\.T+P(c7n:V-O/b7Z4rbreak;}*|&g#w!k/O,F9q//--------------|+P!I8b!_)`r=r*10;for(i=0;in;i++)x_temp=x_min;.e9?9g(|+w-Yf_temp=f_min;S7G&Q3[*`$@}if(f_minf_temp)-I5m!e.H&s%q)y{f_min=f_temp;9A;~;t3lofor(i=0;in;i++)x_min=x_temp;}7Q,Q3f$X+Q:SP9H/K%bf_min=fun(x_min,ys);$kA5},J!Hf_punish=funP(x_min);if(f_ys=1E10)cout\n!!!!!!!!!!!!!!!!!!由于不可知的原因，导致结果发散!!!!!!!!!!!!!!!!!!!endl请在源程序修改合适的精度值或者其他\nendl;cout\n------------外惩罚函数+BFGS变尺度法优化算法C++程序最终结果---------------endl*也可根据实际需要在以上参考值中选择其他结果*endl;cout优化最小值:endl+s1^&s8Y(t/F-|!~.B目标函数f=f_min构造的惩罚函数值=f_punishendl优化程序得到的最优点为:endl|4B1F3c*A3y1|J3jx1=x_min[0]x2=x_min[1]x3=x_min[2]x4=x_min[3]endl约束条件情况:endl约束1=ys[0]+det约束2=ys[1]+det约束3=ys[2]+detendl;}//梯度计算模块//参数：指向目标函数的指针，变量个数，求梯度的点，结果存放处.X+a!^%h`(};I#avoidcomput_grad(double(*pf)(double*x),1a-A*g&~&E7Sintn,/R8l6jS)R,D4}7v.Q8wdouble*point,7O:]2R9P;S)s7{0|8L$^5Ydouble*grad){2T)|8_1@1M'[0Q.odoubleh=1E-5;7l5n#T9l2|inti;double*temp;temp=newdouble[n];for(i=1;i=n;i++)/|-k9[&k9E!V*L{&H8e7g5R+K4K8L0g8}6utemp[i-1]=point[i-1];Q6t!}-I3S&~9M}for(i=1;i=n;i++)/rh4p0ha3B1X9e{)E-Z$y+b#p)Dtemp[i-1]+=0.5*h;grad[i-1]=4*pf(temp)/(3*h);temp[i-1]-=h;/zB7^9L)A.S+s0y*k-O2\grad[i-1]-=4*pf(temp)/(3*h);temp[i-1]+=(3*h/2);grad[i-1]-=(pf(temp)/(6*h));temp[i-1]-=(2*h);grad[i-1]+=(pf(temp)/(6*h));temp[i-1]=point[i-1];}delete[]temp;}//无约束变尺度法BFGS函数声明//参数：pf指向目标函数的指针，n变量个数，min_point接受初始点、存放结果)_%P5[0A*W%{2\5]4v'd//返回0值doubleBFGS(.d9Q$@/u!|6Ddouble(*pf)(double*x),intn,double*min_point)*T,s#c!\5Q&h/G#O{inti,j;intk=0;)q:jT*Yr6_!e;R'Ddoublee=epsilon2;)`/y'w(}6G(S-L4X/Ndoubleg_norm;//梯度的模double*g0=newdouble[n];//梯度double*g1=newdouble[n];double*dg=newdouble[n];double*p=newdouble[n];//搜索方向=-gdoublet;//一维搜索步长8F%Y)^1[#K7Dc'_Sdouble*x0=newdouble[n];4Y;BU#qc's+o%A.Edouble*x1=newdouble[n];double*dx=newdouble[n];double**H=newdouble*[n];for(i=0;in;i++)H=newdouble[n];double*gH=newdouble[n];//计算H变化量的中间变量---double*Hg=newdouble[n];doublenum1;C9}1L)c,F0_6v3s2{doublenum2;//-------------|6pm1q/n#R1T$j+A.j4`;?:s9e//初始化------o#u!Ge0K*\$?*s7qfor(i=0;in;i++)&A${)V!r*Z4y/Ffor(j=0;jn;j++);@Q.D;r;{+j{9M6U,u,j6z-L'[1Hif(i==j)H[j]=1.0;//H0=I2J5x/b'|:H0A$I9{5J%GelseH[j]=0.0;%F-_5h9h2K7D'u/U.`}/t;G*O#W9e.M'g1dfor(i=0;in;i++)-]3F:r.I.X9X,S4I1|#Dx0=min_point;:h6c(k0?3j)Q&qcomput_grad(pf,n,x0,g0);!_2l3?)z/\)Bfor(i=0;in;i++)p=-g0;//--------初始化结束'N.X%U7I-r-x5wdo+S(kr-~&\&?I3J1Q{9V0D-L0be/Bt=line_search(pf,n,x0,p);for(i=0;in;i++)x1=x0+t*p;'{6}.r%k1n3z8f*ocomput_grad(pf,n,x1,g1);//调试程序时参考输出8Y4T7I)J4z5]9b//coutk