惩罚函数法Word下载.docx

1、实验步骤：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.13

2、8664167迭代轮数 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.664902

3、3, 0.5540086, f（ 7）= 0.3111990 0.000003037*外点惩罚函数法优化最优点及目标函数值为： x（ *）= 1.6649023, 0.5540086, f（ *）= 0.3111990算法程序实现#include stdio.hmath.hstdlib.htime.h#define N 2 /*优化设计维数*/#define EPSIN 0.00001 /*迭代精度*/#define H_QJ 1.0 /*初始区间搜索步长*/#define Y_F 1 /*一维搜索方法选择：1黄金分割法*/ /* 2二次插值法*/ #define MC 10.0 /*惩罚因子

4、递增系数*/#define NG 1 /*不等式约束个数*/#define NH 1 /*等式约束个数*/FILE *fp;double mr=1.0; /*外点惩罚因子*/char outname50=外点惩罚函数法计算结果.txt; /*计算结果输出文件*/*给出初始点坐标*/void csd_x（double x0） int i; for（i=0;iN;i+） /*初始点为坐标原点的情况*/ x0i=2.0; return;/*目标函数*/double hanshu1（double x） double f; f=（x0-2）*（x0-2）+（x1-1）*（x1-1）; return f;

5、/*不等式约束方程*/void strain（double x,double g,double h） double eb; eb=EPSIN; g0=-0.25*x0*x0-x1*x1+1.0-eb*10.0;/ g1=1.0-（x0+x1）-eb*10.0; h0=0.0; /*至少有一个为0的等式约束*/*以上为修改部分*/*惩罚函数*/double hanshu（double x） double f,f1,gNG,hNH; f=hanshu1（x）; strain（x,g,h）; f1=0.0;NG;i+） if（gi0.0） f1+=gi*gi;NH; f1+=hi*hi; f+=mr

6、*f1;/*计算f（xk+as）*/double xkadd（double x,double d,double a） double x1N; x1i=xi+a*di; return hanshu（x1）;/*输出选定的一维迭代方法*/void ywddf（int yw） switch（yw） case 1: fprintf（fp,黄金分割法+nn）; break; case 2:二次插值法+nn /*输出当前迭代点坐标及目标函数值*/double xfout（double x,int m） int j; f=hanshu（x）; fprintf（fp, x（%3d）=,m）; for（j=0;

7、jf2） 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 ab0=a1; ab1=a3; else ab0=a3; ab1=a1;/*黄金分割法*/void goldcut（double x,double d,double h,double ebsin） double a1,a2,f1,f2,a,b,ab2; csssqj（x,d,h,ab）; a=ab0; b=ab1; a1=b-0.618*（b-a）; f1=xkadd（x,d,a1）; a2=a+

8、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-aebsin）; xi+=（b+a）*di/2;/*二次插值法*/int rccz（double x,double d,double h,double ebsin） double a1,a2,a3,a4,f1,f2,f3,f4,c1,c2,ab2; int i,p=0,k=0; a1=ab0; a3=ab1;

9、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（a4f4） 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=a

10、2; xi+=f1*di; return p;/*鲍威尔*/void baowr（double x,double h,double ebsin,int yw） double dNN,x0N,x1N,fN+2,df,df1; int i,j,k=1,m; fN=hanshu（x）; for（j=0; dij=0.0; dii=1.0; x0i=xi; f0=fN; switch（yw） case 1: for（i=0; goldcut（x,di,h,ebsin）; fi+1=hanshu（x）; break; case 2: j=rccz（x,di,h,ebsin）; for（i=0; x1i=

11、2*xi-x0i; fN+1=hanshu（x1）; df=0.0; m=0; df1=fi-fi+1; if（df1df） df=df1; m=i; df1=（f0-2*fN+fN+1）*（f0-fN-df）*（f0-fN-df）; df1-=0.5*df*（f0-fN+1）*（f0-fN+1）; if（fN+1=f0|df1 if（fNfN+1） xi=x1i; for（i=m; for（j=0; dij=di+1j; for（i=0; dN-1i=xi-x0i; switch（yw） case 1: goldcut（x,dN-1,h,ebsin）; break; case 2: j=rc

12、cz（x,dN-1,h,ebsin）; fN=hanshu（x）; df+=（xi-x0i）*（xi-x0i）; df=sqrt（df）; k+; while（df/*外点惩罚函数法*/void sumt（double x,double h,double ebsin,int yw） int i,k; double x0N,fact,f,mc; mc=MC;*外点惩罚函数法计算结果*nn鲍威尔法*n ywddf（yw）; /*输出一维迭代方法*/r0=%5.2lf*n,mr）;* 递增系数：c0=%5.2lf*nn,mc）;n f=xfout（x,0）; k=1;迭代轮数 k=%3d n,k）;

13、 baowr（x,h,ebsin,yw）; f=xfout（x,k）; fact=0.0; fact+=（xi-x0i）*（xi-x0i）; fact=sqrt（fact）;%15.9lfnn,fact）; mr*=mc; while（fact*n x（ *）=,xi）;%15.7lf,xN-1）; f（ *）=%15.7lfn,f）;%15.9lfnmain（） double x0N,h,ebsin; int yw; fp=fopen（outname,w csd_x（x0）; h=H_QJ; ebsin=EPSIN; yw=Y_F; sumt（x0,h,ebsin,yw）; /*调用鲍威尔法*/ fclose（fp）; return 0;