变尺度法BFGS算法的C 源码
https://www.wendangku.net/doc/583476211.html,/viewthread.php?tid=73
这个是我那时候交的作业（很粗了。。）。嗯。。。应该标注很清楚了。不过如果你用matlab 编写的话，很多大段的代码，基本上一句话就代替了。
//变尺度法BFGS算法的C++源码
#include "iostream.h"
#include "math.h"
#define A 8
#define B 10
#define C 7
#define D 0
#define START {5,4,5,5} //初始迭代点
const int n=4; //变量个数
const int nc=3; //不等式约束个数(也可根据题目需要增加等式约束个数修改约束构造函数)
const double epsilon1=1E-6; //一维搜索精度
const double epsilon2=1E-4; //变尺度法收敛精度
const double epsilon=1E-4; //惩罚函数法收敛精度
double r=1; //惩罚因子
double br=1; //变尺度法强行退出循环因子(可去除)
const double zr=10; //递增系数
double det=0.0014; //外惩罚函数法约束函数余量
//(此det值仅适合当前A B C D，修改请将其恢复为0.0001或其他建议值)
//计算梯度
void comput_grad(double (*pf)(double *x), int n, double *point, double *grad);
//无约束BFGS变尺度法（目标函数，维数，最小点）返回0值
double BFGS(double (*pf)(double *x), int n, double *min_point);
//构造混合外点惩罚函数返回惩罚函数值（构造点）
double funP(double *p);
//返回目标函数值,并返回约束函数值存入ys （计算点，约束函数值存放处）
double fun(double *x , double *ys );
//一维搜索函数返回最优步长（目标函数，维数，初始点，方向）
double line_search( double (*pf)(double *x),int n,double *start,double *direction);
//配给与一维搜索函数返回下一点的的惩罚函数值（计算点，计算方向，步长）
double funPt(double *x,double *direction,double t);
void main()
{
int i,j;
double c=0; //表示外惩罚函数循环次数
double f_min,f_temp,f_punish;
double *ys=new double[nc];
double x_x=0,f_ys=0;
double x_min[n]=START;
double x_temp[n]=START;
f_temp=fun(x_temp,ys);
for(;;)
{
c++;
BFGS(funP,n,x_min);
f_min=fun(x_min,ys);
for(j=0;j x_x=x_x+pow((x_min[j]-x_temp[j]),2); for(i=0;i f_ys+=fabs(ys); x_x=sqrt(x_x); //参考值显示 cout<<"\n"<<"外惩罚函数循环迭代第 "< <<"目标函数: f_min="< <<"自变量: x1="< x_x="< break;} if((fabs(f_min)<=1)&&(fabs(f_min-f_temp)<=epsilon)) {cout<<"因为迭代xk与xk+1的函数值差距f_f达到收敛精度 "< break;} else if((fabs(f_min)>=1)&&(fabs((f_min-f_temp)/f_min)<=epsilon)) {cout<<"\n 因为迭代xk与xk+1的函数值差距det_f 达到收敛精度"< break;} //--------------| r=r*10; for(i=0;i x_temp=x_min; f_temp=f_min; } if(f_min>f_temp) {f_min=f_temp; for(i=0;i f_min=fun(x_min,ys); f_punish=funP(x_min); if(f_ys>=1E10) cout<<"\n!!!!!!!!!!!!!!!!!!由于不可知的原因，导致结果发散!!!!!!!!!!!!!!!!!!!"< <<"请在源程序修改合适的精度值或者其他\n"< cout<<"\n"<<" ------------外惩罚函数+BFGS变尺度法优化算法C++程序最终结果---------------"< <<" *也可根据实际需要在以上参考值中选择其他结果*"< cout<<"优化最小值:"< <<" 目标函数f="< <<" x1="< <<" 约束1="< }while(g_norm>e); for(i=0;i delete[] g0; delete[] g1; delete[] dg; delete[] p; delete[] x0; delete[] x1; delete[] dx; for (i=0; i delete []H; delete[] gH; delete[] Hg; return 0; } ///////////// double fun(double *x , double *ys )//返回目标函数值 { double f; ys[0]=-x[0]*x[0]-x[1]*x[1]-x[2]*x[2]-x[3]*x[3]-x[0]+x[1]-x[2]+x[3]+A-det; ys[1]=-x[0]*x[0]-2*x[1]*x[1]-x[2]*x[2]-2*x[3]*x[3]+x[0]+x[3]+B-det; ys[2]=-2*x[0]*x[0]-x[1]*x[1]-x[2]*x[2]-2*x[0]+x[1]+x[3]+C-det; f=x[0]*x[0]+x[1]*x[1]+x[2]*x[2]*2+x[3]*x[3]-5*x[0]-5*x[1]-21*x[2]+7*x[3]+D; return f; } double funP(double *p) //构造混合外点惩罚函数返回惩罚函数值{ int i; double f; double *ys=new double[nc]; f=fun(p,ys); if(nc!=0) { for(i=0;i if(ys<0) f+=r*pow(ys[0],2); } return f; } double line_search(double (*pf)(double *x),int n,double *start,double *direction) //参数：指向目标函数的指针，变量个数，出发点，搜索方向 //返回：最优步长 { double a,b,t1=1,t2,t3,tt,h=0.1,f1,f2,f3,j,c=epsilon1; t2=t1+h; f1=funPt(start,direction,t1);f2=funPt(start,direction,t2); //进退法 if(f1>f2) //向前搜索 { t3=t2+h;f3=funPt(start,direction,t3); while(1) if(f2<=f3) {a=t1,b=t3;break;} else{h=2*h;t1=t2;t2=t3;t3=t2+h; f2=funPt(start,direction,t2);f3=funPt(start,direction,t3);} } else { h=-1*h; j=t1;t1=t2;t2=j; j=f1;f1=f2;f2=j; t3=t2+h;f3=funPt(start,direction,t3); while(1) if(f2<=f3){a=t3;b=t1;break;} else{h=2*h;t1=t2;t2=t3;t3=t2+h; f2=funPt(start,direction,t2);f3=funPt(start,direction,t3);} } t1=a+0.382*(b-a); t2=a+0.618*(b-a); f1=funPt(start,direction,t1);f2=funPt(start,direction,t2); while(b-a>c) //黄金分割法， if(f1>f2) {a=t1;t1=t2;t2=a+0.618*(b-a); f1=f2;f2=funPt(start,direction,t2);} else { b=t2;t2=t1;t1=a+0.382*(b-a); f2=f1;f1=funPt(start,direction,t1);} tt=(t1+t2)/2; return tt; } double funPt(double *x,double *direction,double t) //返回带步长为参数的惩罚函数值（用于一维搜索步长时） { double *x1=new double[n]; int i; for(i=0;i x1=x+t*direction; return funP(x1); }