共轭方向法搜索优化c语言源程序

// gong efangxiangfa.cpp : Defines the entry point for the console application.
//

#include
#include
#include

double det=0.0000001; //计算精度
double ak=0.003; //搜索步长


/*目标函数(n维)
入口参数：
x ：n维数组，自变量
返回值 ：函数值 */

double f(double *x)
{
double y;
//y=(x[0]-2)*(x[0]-2)+(x[1]-3)*(x[1]-3)
//y=(x[0]-2)*(x[0]-2)+(x[1]-3)*(x[1]-3)+x[2]*x[2];
//y=(x[0]-2)*(x[0]-2)+(x[1]-3)*(x[1]-3)+x[2]*x[2]+x[3]*x[3];
y=(x[0]-x[1])*(x[0]-x[1])+(x[0]-1)*(x[0]-1)+x[2]*x[2]+x[3]*x[3];
return y;
}

/* 计算X(k+1)=X(k)+a*S(k)
入口参数：
x_k : n维数组，自变量
s_k : n维数组，方向向量
a : 方向步长
n : 维数
出口参数：
x_k1 : 转换后的n维数组*/

void fun(double *x_k1,double *x_k,double *s_k,double a,int n)
{
int i;
for(i=0;ix_k1[i]=x_k[i]+a*s_k[i];
}

/*进退法求谷区间*/
/*入口参数：
n ：优化模型维数
x_k ：n维数组，初始点坐标
s_k ：n维数组，搜索方向向量 */
/*出口参数：
a1 ：搜索区间步长下限
a3 ：搜索区间步长上限*/

void Get_Search_area(int n,double *x_k,double *s_k,double *a1,double *a3)
{
double a=ak,s,f0,a2=0;
double f1=0,f2=0,f3=0;
double *x_k1;
x_k1=(double*)malloc(n*sizeof(double));
f2=f(x_k);
while(1)
{
*a1=a2+a;
fun(x_k1,x_k,s_k,*a1,n);
f1=f(x_k1);
if(f1>f2)
{
if(a==ak)
{
*a3=*a1;
f3=f1;
a=-a;
continue;
}
else if(*a1>*a3)
{
s=*a1;
*a1=*a3;
*a3=s;
f0=f1;
f1=f3;
f3=f0;
break;
}
else break;
}
else
{
*a3=a2;
a2=*a1;
f3=f2;
f2=f1;
a=2*a;
continue;
}
}
}

/*黄金分割法求极小值*/
/*入口参数：
n :优化模型维数
a1 :搜索区间步长下限
a4 :索区间步长上限
s_k :n维数组，搜索方向向量 */
/*出口参数：
x_k :n维数组，极小点坐标 */

void Gold_division(int n,double a1,double a4,double *x_k,double *s_k,double *a_star)
{
double a2,a3,f2,f3;
double *x_k2,*x_k3; int i;
x_k2=(double*)malloc(n*sizeof(double));
x_k3=(double*)malloc(n*sizeof(double));

for(i=0;i{
x_k2[i]=0;
x_k3[i]=0;
}
a2=a1+0.382*(a4-a1);
fun(x_k2,x_k,s_k,a2,n);
f2=f(x_k2);
a3=a1+0.618*(a4-a1);
fun(x_k3,x_k,s_k,a3,n);
f3=f(x_k3);
do
{
if(f2<=f3)
{
a4=a3;
a3=a2;
f3=f2;
a2=a1+0.382*(a4-a1);
fun(x_k2,x_k,s_k,a2,n);
f2=f(x_k2);
}
else
{
a1=a2;
a2=a3;
f2=f3;
a3=a1+0.618*(a4-a1);
fun(x_k3,x_k,s_k,a3,n);
f3=f(x_k3);
}
} while(fabs(a4-a1)>det);
*a_star=(a1+a4)*0.5;
fun(x_k,x_k,s_k,*a_star,n);

free(x_k2);
free(x_k3);
}

/*共轭方向法
/*入口参数：
n :维数
x_k :n维数组，坐标 */
/*出口参数：
f0 :函数极小值 */

void Congjugated_Direction(int n,double *x_k,double *f0)
{
double sum,*f1;
double *a_star,*f_star=0,*x0,a_s=0;
double **s_k,*d,d0=0,*sk,*xk1;
double a1=0,a3=0,f11,f22,f33;
double *a_1,*a_3;
int i,j,m;
a_star=&a_s;
a_1=&a1;
a_3=&a3;
d=(double*)malloc(n*sizeof(double));
x0=(double*)malloc(n*sizeof(double));
f1=(double*)malloc((n+1)*sizeof(double));
xk1=(double*)malloc(n*sizeof(double));
s_k=(double**)malloc(n*sizeof(double*));
sk=(double*)malloc(n*sizeof(double));
for(i=0;i{
s_k[i]=(double*)malloc(n* sizeof(double)); //初始化方向向量矩阵
for(j=0;jif(i==j)
s_k[i][j]=1;
else
s_k[i][j]=0;
}
while(1)
{
f11=f(x_k);
for(i=0;ix0[i]=x_k[i];
f1[0]=f(x_k);
for(i=0;i{
Get_Search_area(n,x_k,s_k[i],a_1,a_3);
Gold_division(n,*a_1,*a_3,x_k,s_k[i],a_star);
f1[i+1]=f(x_k);
d[i]=f1[i]-f1[i+1];
}
f22=f(x_k);
d0=d[0];m=0;printf("a");
for(i=1;iif(d0{
d0=d[i];
m=i;
}
for(i=0;i{
sk[i]=x_k[i]-x0[i];
xk1[i]=2*x_k[i]-x0[i];
}
f33=f(xk1);
if((f33{
Get_Search_area(n,x_k,sk,a_1,a_3);
Gold_division(n,*a_1,*a_3,x_k,sk,a_star);
fun(x_k,x_k,sk,*a_star,n);
for(i=m;ifor(j=0;js_k[m][i]=s_k[m+1][i];
for(j=0;js_k[n-1][j]=sk[j];
}
else
{
if(f22>f33)
{
for(i=0;ix_k[i]=xk1[i];
continue;
}
}
sum=0; //循环终止条件
for(i=0;isum+=(x0[i]-x_k[i])*(x0[i]-x_k[i]);
if(sqrt(sum)<=det)
break;
if(fabs((f(x_k)-f(x0))/f(x0))break;
}
*f0=f(x_k);
free(x0); //释放动态数组
for(i=0;i{
free(s_k[i]);
}
free(s_k);
free(sk);
free(f1);
free(d);
}

/* 主函数*/

int main()
{
double *x_k,*f;
double f0=0;
f=&f0;
int n,i;

puts("原函数为f=(x[0]-x[1])*(x[0]-x[1])+(x[0]-1)*(x[0]-1)+x[2]*x[2]+x[3]*x[3]\n");
puts("请输入变量维数n:\n");
scanf("%d",&n);

x_k=(double*)malloc(n*sizeof(double));

printf("输入起始坐标x_k:\n");
for(i=0;iscanf("%lf",&x_k[i]);

Congjugated_Direction(n,x_k,f);

puts("函数极小点坐标:\n");
for(i=0;iprintf("%lf;",x_k[i]);
printf("\n");

puts("函数最优值:\n");
printf("%lf\n",*f);

free(x_k);
return 0;
}