3.cpp

变尺度法求最优化问题的举例

#include<iostream.h> 
#include<math.h> 
#include<iomanip.h> 
const int n=2; 
static double Rk=10;//Rk是惩罚因子 
static int p=1; 
double fun_original(double x[n]);//既是主函数，也是F(X,Rk)之一 
double a_e_f(double x[n]); 
double fun(double x[n]);//让单纯形法求解此函数最小值 
void g_fun(double x[n],double g0[n]);//牛顿法子函数开始,求出与Rk对应的梯度向量. 
void Hesse_fun(double x[n],double Hesse[n][n]);//求出与Rk对应的Hesse矩阵 
int H(double g0[n],double c); 
void inv(double a[n][n],double c[n][n]);//求Hesse矩阵的逆矩阵 
void xiang_cheng(double a[n][n],double b[n],double c[n]);//n*n矩阵乘以n*1矩阵 
void newton(double x[n]);//牛顿法 
void main() 
{ 
int k=1; 
int i; 
double c_r=0.5;//c是惩罚因子的缩小系数 
double cc=0.000001;//cc是内点法求最终最优解的终止限 
double x[n]={2,2};//给定的初始点在某个增广目标函数的求解域之中 
do 
{ 
newton(x); 
if(fabs(a_e_f(x))<cc)//小心它的值为负，要加上绝对值 
break; 
Rk=c_r*Rk; 
}while(1); 
cout<<"用内点法求解："<<endl; 
cout<<"minf(X)=(X1+1)^3/3+X2"<<endl; 
cout<<" g1(X)=X1-1>=0"<<endl<<"s.t."<<endl;; 
cout<<" g2(X)=X2>=0"<<endl; 
cout<<"结果是："<<endl; 
cout<<"Xk="<<endl; 
for(i=0;i<n;i++) 
cout<<x[i]<<" "; 
cout<<endl; 
cout<<"f(Xk)="<<fun_original(x)<<endl; 
//cout<<"Rk="<<Rk<<endl; 
} 
double fun_original(double x[n])//既是主函数，也是F(X,Rk)之一 
{ 
double y; 
y=pow(x[0]+1,3)/3+x[1]; 
return y; 
} 
double a_e_f(double x[n]) 
{ 
double y; 
y=Rk*(1/(x[0]-1)+1/x[1]); 
return y; 
} 
double fun(double x[n])//让单纯形法求解此函数最小值 
{ 
double y; 
y=fun_original(x)+a_e_f(x); 
return y; 
} 
void g_fun(double x[n],double g0[n])//牛顿法子函数开始,求出与Rk对应的梯度向量. 
{ 
g0[0]=pow(x[0]+1,2)-Rk/pow(x[1]-1,2); 
g0[1]=1-Rk/pow(x[1],2); 
} 
void Hesse_fun(double x[n],double Hesse[n][n])//求出与Rk对应的Hesse矩阵 
{ 
Hesse[0][0]=2*(x[0]+1)+2*Rk/pow(x[1]-1,3); 
Hesse[0][1]=Hesse[1][0]=0; 
Hesse[1][1]=2*Rk/pow(x[1],3); 
} 
int H(double g0[n],double c) 
{ 
double s=0; 
for(int i=0;i<n;i++) 
{ 
s+=pow(g0[i],2); 
} 
if(sqrt(s)<c) 
return 1; 
else 
return 0; 
} 
void inv(double a[n][n],double c[n][n])//求Hesse矩阵的逆矩阵 
{ 
double m[n];//辅助乘数 
double T;//存储换行临时变量 
double temp=0; 
double t;//最大列主元 
int tap1=0,tap2=0;//最大列主元下标 
for(int i=0;i<n;i++) 
{ 
for(int j=0;j<n;j++) 
{ 
if(i==j) 
{ 
c[i][j]=1; 
} 
else 
{ 
c[i][j]=0; 
} 
} 
} 
for(int s=0;s<n;s++) 
{ 
t=a[s][s];//赋初值 
for(int p=s;p<n;p++)//选取最大列主元 
{ 
if(fabs(a[p][s])>t) 
{ 
t=a[p][s]; 
tap1=p; 
tap2=s; 
} 
} 
if(t==0) 
{cout<<"列主元为零，失败！"<<endl;break;} 
if(tap1!=tap2)//换行 
{ 
for(int i=0;i<n;i++)////////////////////////////////s 
{ 
T=a[tap1][i]; 
a[tap1][i]=a[s][i]; 
a[s][i]=T; 
T=c[tap1][i];//逆矩阵换行 
c[tap1][i]=c[s][i]; 
c[s][i]=T; 
} 
} 
tap1=tap2=0;//置零 
for(int j=0;j<n;j++)//单位化,(注意：逆矩阵每列都单位化，而原矩阵则不受影响) 
{ 
a[s][j]=a[s][j]/t; 
c[s][j]=c[s][j]/t; 
} 
for(int i=0;i<n;i++)//消元 
{ 
if(i!=s) 
{ 
m[i]=a[i][s]/a[s][s]; 
for(int j=0;j<n;j++)//注意：逆矩阵每列都消元，而原矩阵则不受影响 
{ 
a[i][j]=a[i][j]-m[i]*a[s][j]; 
c[i][j]=c[i][j]-m[i]*c[s][j]; 
} 
} 
} 
} 
} 
void xiang_cheng(double a[n][n],double b[n],double c[n])//n*n矩阵乘以n*1矩阵 
{ 
int i,j,k; 
if(n==n) 
{ 
//cout<<"两矩阵阶数相符可以相乘！"<<endl; 
for(i=0;i<n;i++) 
{ 
c[i]=0; 
} 
for(k=0;k<n;k++)//c矩阵的第k行 
{ 
for(j=0;j<n;j++) 
{ 
c[k]+=a[k][j]*b[j]; 
} 
} 
} 
else 
cout<<"两个矩阵阶数不符，不能相乘"<<endl; 
} 
void newton(double x[n]) 
{ 
int i; 
double g0[n];//梯度的值 
double c_H=0.000001;//H准则限值 
double Hesse[n][n];//Hesse矩阵 
double inv_Hesse[n][n];//Hesse矩阵的逆 
double temp[n];//保存矩阵相乘结果 

g_fun(x,g0);//返回梯度的初值 
do 
{ 
Hesse_fun(x,Hesse);//返回Hesse矩阵的值？？一个二次函数的Hesse矩阵是变的还是不变 
inv_Hesse[0][1]=inv_Hesse[1][0]=0; 
inv_Hesse[0][0]=1/Hesse[0][0]; 
inv_Hesse[1][1]=1/Hesse[1][1]; 
xiang_cheng(inv_Hesse,g0,temp); 
for(i=0;i<n;i++) 
{ 
x[i]=x[i]-temp[i]; 
} 
g_fun(x,g0); 
}while(H(g0,c_H)==0); 
//cout<<"使用次"<<p<<"牛顿法"<<endl; 
p++; 
}