//多维无约束优化软件设计#include

#include

#include

double det=1e-5; //计算精度double det1=1e-3; //梯度判断精度double ak=3e-3; //搜索步长double dx=1e-4; //梯度计算步长

/*目标函数(n维)

入口参数：

x ：n维数组，自变量

返回值 ：函数值*/

double f(double *x)

{

double y;

y=(x[0]-2)*(x[0]-2)*(x[0]-2)*(x[0]-2)+(x[0]-2*x[1])*(x[0]-2*x[1]);

//y=(x[0]-2)^4+(x[0]-2*x[1])^2; 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;i

x_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(sizeof(double)*n);

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;

}

}

free(x_k1);

}

/*黄金分割法求极小值*/

/*入口参数：

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 a2,a3,f2,f3,a_star;

double *x_k2,*x_k3; int i;

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

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

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 ：优化模型维数

dx ：梯度步长

x_k ：n维数组，初始点坐标

出口参数：

s_k ：梯度的负方向向量*/

void tidu(int n,double dx,double x[],double s_k[])

{

int i;

double *x0;

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

for(i=0;i

x0[i]=x[i];

for(i=0;i

{

x0[i]=x[i]+dx;

s_k[i]=(f(x0)-f(x))/dx*(-1);

x0[i]=x0[i]-dx;

}

free(x0);

}

/*计算梯度范数

入口参数：

n ：优化模型维数

s_k ：梯度的负方向向量

返回值：梯度范数*/

double f_fanshu(int n,double s_k[])

{

int i;

double sum=0;

for(i=0;i

sum+=s_k[i]*s_k[i];

return (sqrt(sum));

}

/*梯度法求最优函数值

/*入口参数：

n :优化模型维数

x_k :n维数组，极小点坐标*/

/*出口参数：

f0 :函数极小值*/

void fun_tidu(int n,double *x_k,double *f0)

{

double *s_k;

double a1=0,a3=0;

double *a_1,*a_3;

a_1=&a1;

a_3=&a3;

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

while(1)

{ tidu(n,dx,x_k,s_k);

if(f_fanshu(n,s_k)>det1)

{

Get_Search_area(n,x_k,s_k,a_1,a_3);

Gold_division(n,*a_1,*a_3,x_k,s_k);

}

else

break;

}

*f0=f(x_k);

free(s_k);

}

/*主函数*/

int main()

{

int n,i;

double *x_k,*f;

double f0=0;

f=&f0;

puts("原函数为f=(x[0]-2)*(x[0]-2)*(x[0]-2)*(x[0]-2)+(x[0]-2*x[1])*(x[0]-2*x[1]);\n");

puts("请输入变量维数n:\n");

scanf("%d",&n);

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

printf("输入起始坐标x_k:\n");

for(i=0;i

{

scanf("%lf",&x_k[i]);

fun_tidu(n,x_k,f);

}

puts("函数极小点坐标:\n");

for(i=0;i

printf("%lf;",x_k[i]);

printf("\n");

puts("函数最优值:\n");

printf("%lf\n",*f);

free(x_k);

return 0;

}