惩罚函数法Word下载.docx
《惩罚函数法Word下载.docx》由会员分享，可在线阅读，更多相关《惩罚函数法Word下载.docx（16页珍藏版）》请在冰豆网上搜索。
实验步骤：
1，画流程图，编写程序；
2，将目标函数代入；
3，编译运行，将结果保存。
结
果
分
析
**********外点惩罚函数法计算结果***********
**********无约束优化方法：
鲍威尔法*********
++++++一维搜索方法：
黄金分割法++++++
**********初始惩罚因子：
r0=1.00**********
**********递增系数：
c0=10.00*************
初始坐标：
x（0）=[2.0000000,2.0000000],f（0）=17.0008000
迭代轮数k=1
x
（1）=[1.7819615,0.6713951],f
（1）=0.2154089
迭代精度：
1.346377264
迭代轮数k=2
x
（2）=[1.6850587,0.5722104],f
（2）=0.2961649
0.138664167
迭代轮数k=3
x（3）=[1.6671097,0.5559496],f（3）=0.3095916
0.024219441
迭代轮数k=4
x（4）=[1.6651347,0.5541996],f（4）=0.3110373
0.002638740
迭代轮数k=5
x（5）=[1.6649374,0.5540160],f（5）=0.3111830
0.000269539
迭代轮数k=6
x（6）=[1.6649051,0.5540098],f（6）=0.3111975
0.000032957
迭代轮数k=7
x（7）=[1.6649023,0.5540086],f（7）=0.3111990
0.000003037
*********************
外点惩罚函数法优化最优点及目标函数值为：
x（*）=[1.6649023,0.5540086],f（*）=0.3111990
算法程序实现
#include<
string.h>
stdio.h>
math.h>
stdlib.h>
time.h>
#defineN2/*优化设计维数*/
#defineEPSIN0.00001/*迭代精度*/
#defineH_QJ1.0/*初始区间搜索步长*/
#defineY_F1/*一维搜索方法选择：
1——黄金分割法*/
/*2——二次插值法*/
#defineMC10.0/*惩罚因子递增系数*/
#defineNG1/*不等式约束个数*/
#defineNH1/*等式约束个数*/
FILE*fp;
doublemr=1.0;
/*外点惩罚因子*/
charoutname[50]="
外点惩罚函数法计算结果.txt"
;
/*计算结果输出文件*/
/*给出初始点坐标*/
voidcsd_x（doublex0[]）
{
inti;
for（i=0;
i<
N;
i++）/*初始点为坐标原点的情况*/
x0[i]=2.0;
return;
}
/*目标函数*/
doublehanshu1（doublex[]）
doublef;
f=（x[0]-2）*（x[0]-2）+（x[1]-1）*（x[1]-1）;
returnf;
/*不等式约束方程*/
voidstrain（doublex[],doubleg[],doubleh[]）
doubleeb;
eb=EPSIN;
g[0]=-0.25*x[0]*x[0]-x[1]*x[1]+1.0-eb*10.0;
//g[1]=1.0-（x[0]+x[1]）-eb*10.0;
h[0]=0.0;
/*至少有一个为0的等式约束*/
/***********以上为修改部分***********/
/*惩罚函数*/
doublehanshu（doublex[]）
doublef,f1,g[NG],h[NH];
f=hanshu1（x）;
strain（x,g,h）;
f1=0.0;
NG;
i++）
if（g[i]<
0.0）
f1+=g[i]*g[i];
NH;
f1+=h[i]*h[i];
f+=mr*f1;
/*计算f（xk+as）*/
doublexkadd（doublex[],doubled[],doublea）
doublex1[N];
x1[i]=x[i]+a*d[i];
returnhanshu（x1）;
/*输出选定的一维迭代方法*/
voidywddf（intyw）
switch（yw）
{
case1:
fprintf（fp,"
黄金分割法++++++\n\n"
）;
break;
case2:
二次插值法++++++\n\n"
}
/*输出当前迭代点坐标及目标函数值*/
doublexfout（doublex[],intm）
intj;
f=hanshu（x）;
fprintf（fp,"
x（%3d）=["
m）;
for（j=0;
j<
N-1;
j++）
fprintf（fp,"
%15.7lf,"
x[j]）;
%15.7lf],f（%3d）=%15.7lf\n"
x[N-1],m,f）;
/*初始搜索区间的确定*/
voidcsssqj（doublex[],doubled[],doubleh,doubleab[]）
doublea1,a2,a3,f1,f2,f3;
a2=0.0;
a3=a2+h;
f2=xkadd（x,d,a2）;
f3=xkadd（x,d,a3）;
if（f3>
f2）
a2=a3;
a3=0.0;
f1=f2;
f2=f3;
f3=f1;
h=-h;
do
a1=a2;
a3=a2+h;
f3=xkadd（x,d,a3）;
h=2*h;
}while（f3<
f2）;
if（h>
ab[0]=a1;
ab[1]=a3;
else
ab[0]=a3;
ab[1]=a1;
/*黄金分割法*/
voidgoldcut（doublex[],doubled[],doubleh,doubleebsin）
doublea1,a2,f1,f2,a,b,ab[2];
csssqj（x,d,h,ab）;
a=ab[0];
b=ab[1];
a1=b-0.618*（b-a）;
f1=xkadd（x,d,a1）;
a2=a+0.618*（b-a）;
if（f1>
{
a=a1;
a1=a2;
f1=f2;
a2=a+0.618*（b-a）;
f2=xkadd（x,d,a2）;
}
else
b=a2;
a2=a1;
f2=f1;
a1=b-0.618*（b-a）;
f1=xkadd（x,d,a1）;
}while（b-a>
ebsin）;
x[i]+=（b+a）*d[i]/2;
/*二次插值法*/
intrccz（doublex[],doubled[],doubleh,doubleebsin）
doublea1,a2,a3,a4,f1,f2,f3,f4,c1,c2,ab[2];
inti,p=0,k=0;
a1=ab[0];
a3=ab[1];
a2=（a1+a3）/2;
while
（1）
c1=（f3-f1）/（a3-a1）;
c2=（（f2-f1）/（a2-a1）-c1）/（a2-a3）;
if（0==c2）
p=1;
a4=0.5*（a1+a3-c1/c2）;
if（（a4-a1）*（a3-a4）<
=0.0）
p=2;
f4=xkadd（x,d,a4）;
if（1==k&
&
fabs（a4-a2）<
=ebsin）
if（a4<
a2）
if（f2>
f4）
{
f1=f2;
a1=a2;
k=1;
a2=a4;
f2=f4;
}
else
a3=a4;
f3=f4;
f3=f2;
a3=a2;
a1=a4;
f1=f4;
if（0==p）
f1=f4-f2?
a4:
a2;
f1=a2;
x[i]+=f1*d[i];
returnp;
/*鲍威尔*/
voidbaowr（doublex[],doubleh,doubleebsin,intyw）
doubled[N][N],x0[N],x1[N],f[N+2],df,df1;
inti,j,k=1,m;
f[N]=hanshu（x）;
for（j=0;
d[i][j]=0.0;
d[i][i]=1.0;
x0[i]=x[i];
f[0]=f[N];
switch（yw）
case1:
for（i=0;
{
goldcut（x,d[i],h,ebsin）;
f[i+1]=hanshu（x）;
}
break;
case2:
j=rccz（x,d[i],h,ebsin）;
}
for（i=0;
x1[i]=2*x[i]-x0[i];
f[N+1]=hanshu（x1）;
df=0.0;
m=0;
df1=f[i]-f[i+1];
if（df1>
df）
df=df1;
m=i;
df1=（f[0]-2*f[N]+f[N+1]）*（f[0]-f[N]-df）*（f[0]-f[N]-df）;
df1-=0.5*df*（f[0]-f[N+1]）*（f[0]-f[N+1]）;
if（f[N+1]>
=f[0]||df1>
if（f[N]>
f[N+1]）
x[i]=x1[i];
for（i=m;
for（j=0;
d[i][j]=d[i+1][j];
for（i=0;
d[N-1][i]=x[i]-x0[i];
switch（yw）
case1:
goldcut（x,d[N-1],h,ebsin）;
break;
case2:
j=rccz（x,d[N-1],h,ebsin）;
f[N]=hanshu（x）;
df+=（x[i]-x0[i]）*（x[i]-x0[i]）;
df=sqrt（df）;
k++;
}while（df>
/*外点惩罚函数法*/
voidsumt（doublex[],doubleh,doubleebsin,intyw）
inti,k;
doublex0[N],fact,f,mc;
mc=MC;
**********外点惩罚函数法计算结果**********\n\n"
鲍威尔法**********\n"
ywddf（yw）;
/*输出一维迭代方法*/
r0=%5.2lf**********\n"
mr）;
**********递增系数：
c0=%5.2lf**********\n\n"
mc）;
\n"
f=xfout（x,0）;
k=1;
迭代轮数k=%3d\n"
k）;
baowr（x,h,ebsin,yw）;
f=xfout（x,k）;
fact=0.0;
fact+=（x[i]-x0[i]）*（x[i]-x0[i]）;
fact=sqrt（fact）;
%15.9lf\n\n"
fact）;
mr*=mc;
}while（fact>
*********************\n"
x（*）=["
x[i]）;
%15.7lf],"
x[N-1]）;
f（*）=%15.7lf\n"
f）;
"
%15.9lf\n"
main（）
doublex0[N],h,ebsin;
intyw;
fp=fopen（outname,"
w"
csd_x（x0）;
h=H_QJ;
ebsin=EPSIN;
yw=Y_F;
sumt（x0,h,ebsin,yw）;
/*调用鲍威尔法*/
fclose（fp）;
return0;
升级会员 