外点罚函数法

太原理工大学机械学院机测系课程上机实验报告课程名称：机械优化设计班级机测2班日期成绩评定姓名陈明鑫实验室机械馆机房老师签名实验名称利用外点罚函数法求解相关函数的极小值点所用软件实验目的及内容实验目的：1.掌握并能够建立最优化基本类型问题的数学模型。2.掌握最优化方法的基本概念、基本理论和基本方法，奠定最优化的理论基础。3.能够熟练编制和调试最优化方法的程序，奠定解决实际中的优化问题的基础实验内容：理解外点罚函数法并编写相关程序求其极小值点。实验原理步骤、实验结果及分析外点罚函数法计算结果**********外点惩罚函数法计算结果********************无约束优化方法：鲍威尔法**********++++++一维搜索方法：黄金分割法++++++**********初始惩罚因子：r0=1.00********************递增系数：c0=10.00**********初始坐标：x(0)=[0.0000000,0.0000000],f(0)=1.0000200迭代轮数k=1x(1)=[0.5000044,0.0000007],f(1)=0.5000100迭代精度：0.500004358迭代轮数k=2x(2)=[0.9090993,-0.0000014],f(2)=0.9091091迭代精度：0.409094982迭代轮数k=3x(3)=[0.9901088,0.0000007],f(3)=0.9901188迭代精度：0.081009509迭代轮数k=4x(4)=[0.9990110,-0.0000008],f(4)=0.9990210迭代精度：0.008902129迭代轮数k=5x(5)=[0.9999100,0.0000007],f(5)=0.9999200迭代精度：0.000899031迭代轮数k=6x(6)=[1.0000000,-0.0000008],f(6)=1.0000100迭代精度：0.000090004迭代轮数k=7x(7)=[1.0000090,0.0000006],f(7)=1.0000190迭代精度：0.000009112迭代轮数k=8x(8)=[1.0000099,-0.0000007],f(8)=1.0000199迭代精度：0.000001643迭代轮数k=9x(9)=[1.0000100,-0.0000011],f(9)=1.0000200迭代精度：0.000000346*********************外点惩罚函数法优化最优点及目标函数值为：x(*)=[1.0000100,-0.0000011],f(*)=1.0000200迭代精度：0.000000346罚函数法考核题外点：0120125.0..)1()2()(min2122212221xxxxtsxxxf，其初始迭代点为Tx]2,2[0，中止误差51e内点：0120125.0..)1()2()(min2122212221xxxxtsxxxf，其初始迭代点为Tx]2,2[0，中止误差51e**********外点惩罚函数法计算结果********************无约束优化方法：鲍威尔法**********++++++一维搜索方法：黄金分割法++++++**********初始惩罚因子：r0=1.00********************递增系数：c0=10.00**********初始坐标：x(0)=[2.0000000,2.0000000],f(0)=18.0000000迭代轮数k=1x(1)=[1.3416070,0.9778107],f(1)=0.7478686迭代精度：1.215875094迭代轮数k=2x(2)=[0.9289099,0.9293352],f(2)=1.2645789迭代精度：0.415534359迭代轮数k=3x(3)=[0.8351188,0.9136297],f(3)=1.3788208迭代精度：0.095096941迭代轮数k=4x(4)=[0.8241206,0.9116622],f(4)=1.3919793迭代精度：0.011172862迭代轮数k=5x(5)=[0.8230004,0.9114603],f(5)=1.3933162迭代精度：0.001138249迭代轮数k=6x(6)=[0.8228881,0.9114401],f(6)=1.3934501迭代精度：0.000114040迭代轮数k=7x(7)=[0.8228790,0.9114382],f(7)=1.3934679迭代精度：0.000009338迭代轮数k=8x(8)=[0.8228774,0.9114380],f(8)=1.3934925迭代精度：0.000001571迭代轮数k=9x(9)=[0.8228757,0.9114378],f(9)=1.3934650迭代精度：0.000001763迭代轮数k=10x(10)=[0.8228757,0.9114378],f(10)=1.3934650迭代精度：0.000000011*********************外点惩罚函数法优化最优点及目标函数值为：x(*)=[0.8228757,0.9114378],f(*)=1.3934650迭代精度：0.000000011算法程序实现//*鲍威尔法*/#includestring.h#includestdio.h#includemath.h#includestdlib.h#includetime.h#defineN2/*优化设计维数*/#defineEPSIN0.000001/*迭代精度*/#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;iN;i++)/*初始点为坐标原点的情况*/x0[i]=0.0;return;}/*目标函数*/doublehanshu1(doublex[]){doublef;f=x[0]*x[0]+x[1]*x[1];returnf;}/*不等式约束方程*/voidstrain(doublex[],doubleg[],doubleh[]){doubleeb;eb=EPSIN;g[0]=x[0]-1.0-eb*10.0;//g[1]=1.0-(x[0]+x[1])-eb*10.0;h[0]=0.0;/*至少有一个为0的等式约束*/return;}/***********以上为修改部分***********//*惩罚函数*/doublehanshu(doublex[]){inti;doublef,f1,g[NG],h[NH];f=hanshu1(x);strain(x,g,h);f1=0.0;for(i=0;iNG;i++)if(g[i]0.0)f1+=g[i]*g[i];for(i=0;iNH;i++)f1+=h[i]*h[i];f+=mr*f1;returnf;}/*计算f(xk+as)*/doublexkadd(doublex[],doubled[],doublea){inti;doublex1[N];for(i=0;iN;i++)x1[i]=x[i]+a*d[i];returnhanshu(x1);}/*输出选定的一维迭代方法*/voidywddf(intyw){switch(yw){case1:fprintf(fp,++++++一维搜索方法：黄金分割法++++++\n\n);break;case2:fprintf(fp,++++++一维搜索方法：二次插值法++++++\n\n);break;}return;}/*输出当前迭代点坐标及目标函数值*/doublexfout(doublex[],intm){intj;doublef;f=hanshu(x);fprintf(fp,x(%3d)=[,m);for(j=0;jN-1;j++)fprintf(fp,%15.7lf,,x[j]);fprintf(fp,%15.7lf],f(%3d)=%15.7lf\n,x[N-1],m,f);returnf;}/*初始搜索区间的确定*/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(f3f2){a2=a3;a3=0.0;f1=f2;f2=f3;f3=f1;h=-h;}do{a1=a2;a2=a3;f1=f2;f2=f3;a3=a2+h;f3=xkadd(x,d,a3);h=2*h;}while(f3f2);if(h0.0){ab[0]=a1;ab[1]=a3;}else{ab[0]=a3;ab[1]=a1;}return;}/*黄金分割法*/voidgoldcut(doublex[],doubled[],doubleh,doubleebsin){doublea1,a2,f1,f2,a,b,ab[2];inti;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);f2=xkadd(x,d,a2);do{if(f1f2){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);for(i=0;iN;i++)x[i]+=(b+a)*d[i]/2;return;}/*二次插值法*/intrccz(doublex[],doubled[],doubleh,doubleebsin){doublea1,a2,a3,a4,f1,f2,f3,f4,c1,c2,ab[2];inti,p=0,k=0;csssqj(x,d,h,ab);a1=ab[0];a3=ab[1];a2=(a1+a3)/2;f1=xkadd(x,d,a1);f2=xkadd(x,d,a2);f3=xkadd(x,d,a3);while(1){c1=(f3-f1)/(a3-a1);c2=((f2-f1)/(a2-a1)-c1)/(a2-a3);if(0==c2){p=1;break;}a4=0.5*(a1+a3-c1/c2);if((a4-a1)*(a3-a4)=0.0){p=2;break;}f4=xkadd(x,d,a4);if(1==k&&fabs(a4-a2)=ebsin)break;if(a4a2){if(f2f4){f1=f2;a1=a2;k=1;a2=a4;f2=f4;}else{a3=a4;f3=f4;}}else{if(f2f4){f3=f2;a3=a2;k=1;a2=a4;f2=f4;}else{a1=a4;f1=f4;}}}if(0==p)f1=f4-f2?a4:a2;elsef1=a2;for(i=0;iN;i++)x[i]+=f1*d[i];returnp;}/*鲍威尔*/voidbaowr(doublex[],doubleh,doubleebsin,intyw){doubled[N][N],x0