混合法求极值程序C语言

优化设计作业混合惩罚函数法求极值轮机研李凯2011129026
混合惩罚函数法求极值程序
#include<stdio.h>
#include<math.h>
double dv[2];
double x1[2];
//目标函数
double func(double x1,double x2)
{
return ((x1-2)*(x1-2)+(x2-1)*(x2-1));
}
//约束条件
double hj(double x1,double x2)
{
return (x1-2*x2+1);
}
double gi(double x1,double x2)
{
double g=0.25*x1*x1+x2*x2-1;
if (g<0)
return 0;
else return g;
}
//惩罚函数
double fun(double x1,double x2,double r) {
double f=func(x1,x2);
double h=hj(x1,x2);
double g=gi(x1,x2);
return (f+r/g+(g*g+h*h)/r);
}
//单纯形法求最优解
void F(double x1,double x2,double dv[2],double r)
{
double x[7][2];
double max,min,t,xh[2],xe[2],f[6],fe,fh,fl; int i,l,m,n,k=0;
x[0][0]=x1;
x[0][1]=x2;
x[1][0]=2.5;
x[1][1]=3;
x[2][0]=3;
x[2][1]=2.5;
do
{
for(i=0;i<=2;i++)
f[i]=fun(x[i][0],x[i][1],r);
/*****判断好点、坏点*****/
max=min=f[0];
m=n=0;
for(i=0;i<=2;i++)
{
if(f[i]>max)
{max=f[i]; m=i;}
if(f[i]<min)
{min=f[i]; n=i;}
}
xh[0]=x[m][0];
xh[1]=x[m][1];
fh=fun(xh[0],xh[1],r);
xe[0]=x[n][0];
xe[1]=x[n][1];
fe=fun(xe[0],xe[1],r);
/*****求中心点、反射点*****/
for(i=0;i<2;i++)
{
x[3][i]=0.5*(x[0][i]+x[1][i]+x[2][i]-xh[i] );
x[4][i]=x[3][i]+1*(x[3][i]-xh[i]);
}
f[3]=fun(x[3][0],x[3][1],r);
f[4]=fun(x[4][0],x[4][1],r);
/*****判断除好点、坏点以外的点l*****/
for(l=0;l<=2;l++)
if(l!=m&&l!=n)break;
fl=fun(x[l][0],x[l][1],r);
/*****分类判断*****/
if(f[4]>fe) /*****反射点函数值比好点差*****/
if(f[4]>fl) /*****反射点函数值比除坏点以外的点都差*****/
{ if(f[4]>=fh) /*****反射点函数值比坏点还差，反射失败*****/
{
x[6][0]=x[3][0]+0.5*(xh[0]-x[3][0]);
x[6][1]=x[3][1]+0.5*(xh[1]-x[3][1]);
f[6]=fun(x[6][0],x[6][1],r);
if(f[6]>fh) /*****收缩失败，缩小边长*****/
for(i=0;i<=2;i++)
{
x[i][0]=0.5*(x[i][0]+xe[0]);
x[i][1]=0.5*(x[i][1]+xe[1]);
}
else /*****收缩成功，形成新的单纯矩阵*****/
x[m][0]=x[6][0];
x[m][1]=x[6][1];
}
else /*****反射点仅比最坏点好*****/
{
xh[0]=x[4][0];xh[1]=x[4][1];
x[6][0]=x[3][0]+0.5*(xh[0]-x[3][0]);
x[6][1]=x[3][1]+0.5*(xh[1]-x[3][1]);
f[6]=fun(x[6][0],x[6][1],r);
if(f[6]>fh)
for(i=0;i<=2;i++)
{
x[i][0]=0.5*(x[i][0]+ xe[0]);
x[i][1]=0.5*(x[i][1]+ xe[1]);
}
else
x[m][0]=x[6][0];
x[m][1]=x[6][1];
}
}
else /*****反射点比除坏点以外的点好*****/
{ x[m][0]=x[4][0];
x[m][1]=x[4][1];
}
else /*****反射点比好点还好*****/
{ x[5][0]=x[3][0]+2*(x[4][0]-x[3][0]);
x[5][1]=x[3][1]+2*(x[4][1]-x[3][1]);
f[5]=fun(x[5][0],x[5][1],r);
if(f[5]>fe)
{ x[m][0]=x[4][0]; x[m][1]=x[4][1]; }
else
{ x[m][0]=x[5][0];
x[m][1]=x[5][1]; }
}
for(i=0;i<=2;i++)
{
t=fabs(f[0]-f[3]);
if(fabs(f[i]-f[3])>=t)t=fabs(f[i]-f[3]);
else break;
}
}while(t>1e-4);
dv[0]=x[n][0];
dv[1]=x[n][1];
printf("dv0=%f %f\n",dv[0],dv[1]); }
main()
{
int n=1,c=10,i;
double x[2]={1.5,2};
double e=0.001,f1,f2,r=0.00000001,w1;
do
{
for(i=0;i<2;i++)
x1[i]=x[i];
F(x[0],x[1],dv,r);
printf("dv=%f %f\n",dv[0],dv[1]);
x[0]=dv[0];
x[1]=dv[1];
printf("\n\nr=%f,x1=%f,x2=%f\n",r,x[0], x[1]);
if(fabs(sqrt((x[0]-x1[0])*(x[0]-x1[0])+(x[ 1]-x1[0])*(x[1]-x1[1])))>=1e-3)
{r=c*r;
printf("r=%f ",r);
}
else break;
}while(1);
w1=func(x[0],x[1]);
printf("x1=%f,x2=%f,w1=%f
",x[0],x[1],w1);
}
运行结果：。