回复 2楼 beyondyf
通过进退法确定谷区间和黄金分割法求极小值然后求极小值，我现在将我的代码发给您看看。前段时间因为出了点事，所以拖到了现在，万分抱歉。
程序代码：
// 多维无约束优化软件设计
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
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<n;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<n;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<n;i++)
x0[i]=x[i];
for(i=0;i<n;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<n;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<n;i++)
{
scanf("%lf",&x_k[i]);
fun_tidu(n,x_k,f);
}
puts("函数极小点坐标:\n");
for(i=0;i<n;i++)
printf("%lf;",x_k[i]);
printf("\n");
puts("函数最优值:\n");
printf("%lf\n",*f);
free(x_k);
return 0;
}
2014-02-23 21:49
