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

//变尺度法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;jx_x=x_x+pow((x_min[j]-x_temp[j]),2);
for(i=0;if_ys+=fabs(ys);
x_x=sqrt(x_x);
//参考值显示
cout<<"\n"<<"外惩罚函数循环迭代第 "<<<"目标函数: f_min="<<<"自变量: x1="<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;ix_temp=x_min;
f_temp=f_min;
}
if(f_min>f_temp)
{f_min=f_temp;
for(i=0;if_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;idelete[] g0;
delete[] g1;
delete[] dg;
delete[] p;
delete[] x0;
delete[] x1;
delete[] dx;
for (i=0; idelete []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;iif(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;ix1=x+t*direction;
return funP(x1);

}