无约束优化方法---鲍威尔方法本实验用鲍威尔方法求函数f(x)=(x1-5)2+(x2-6)2 的最优解。
一、简述鲍威尔法的基本原理从任选的初始点x⑴o出发，先按坐标轮换法的搜索方向依次沿e1.e2.e3进行一维搜索，得各自方向的一维极小点x⑴ x⑵ x⑶.连接初始点xo⑴和最末一个一维极小点x3⑴,产生一个新的矢量 S1=x3⑴-xo⑴再沿此方向作一维搜索，得该方向上的一维极小点x⑴.从xo⑴出发知道获得x⑴点的搜索过程称为一环。
S1是该环中产生的一个新方向，称为新生方向。
接着，以第一环迭代的终点x⑴作为第二环迭代的起点xo⑵,即 Xo⑵←x⑴弃去第一环方向组中的第一个方向e1，将第一环新生方向S1补在最后，构成第二环的基本搜索方向组e2,e3,S1，依次沿这些方向求得一维极小点x1⑵,x2⑵,x3⑵.连接Xo⑵与x3⑵,又得第二环的新生方向S2=x3⑵-xo⑵沿S2作一维搜索所得的极小点x⑵即为第二环的最终迭代点二、鲍威尔法的程序#include "stdafx.h" /* 文件包含*/#include <math.h>#include <stdio.h>#include <stdlib.h>#define MAXN 10#define sqr(x) ((x)*(x))double xkk[MAXN],xk[MAXN],sk[MAXN];int N,type,nt,et; //N--变量个数，type=0,1,2,3 nt,et--不等式、等式约束个数double rk;double funt(double *x,double *g,double *h){g[0]=x[0]; g[1]=x[1]-1; g[2]=11-x[0]-x[1];return sqr(x[0]-8)+sqr(x[1]-8);}double F(double *x){double f1,f2,ff,fx,g[MAXN],h[MAXN];int i;fx=funt(x,g,h);f1=f2=0.0;if(type==0 || type==2)for(i=0; i<nt; i++) f1+=(fabs(g[i])>1.0e-15)?1.0/g[i]:1.0e15;else for(i=0; i<nt; i++) f1+=(g[i]>0)?0:g[i]*g[i];for(i=0; i<et; i++) f2+=h[i]*h[i];if(type==0 || type==2)ff=fx+rk*f1+f2/sqrt(rk);else ff=fx+rk*(f1+f2);return ff;}double f(double x){int i;for (i=0; i<N; i++) xkk[i]=xk[i]+x*sk[i];return F(xkk);}void find_ab(double x0,double h0,double *a,double *b) {double h,x1,y1,x2,y2,x3,y3;h=h0;x1=x0; y1=f(x1); x2=x1+h; y2=f(x2);if (y2>=y1) {h=-h0; x3=x1; y3=y1;x1=x2; y1=y2; x2=x3; y2=y3;}for (;;) {h*=2.0; x3=x2+h; y3=f(x3);if (y2<y3) break;x1=x2; y1=y2; x2=x3; y2=y3;}if (h>0) {*a=x1; *b=x3;}else {*a=x3; *b=x1;}}void search_gold(double a,double b,double e,double *x,double *y) {double x1,x2,y1,y2;x1=a+0.382*(b-a); y1=f(x1);x2=a+0.618*(b-a); y2=f(x2);do {if (y1<y2) {b=x2; x2=x1; y2=y1;x1=a+0.382*(b-a); y1=f(x1);} else {a=x1; x1=x2; y1=y2;x2=a+0.618*(b-a); y2=f(x2);}} while ((b-a)>e);*x=0.5*(a+b); *y=f(*x);}double nc_powell(double *x0,double e1,double e2){int i,j,k=1,m;double a,b,ax,ay,d;doubless[MAXN][MAXN],s1[MAXN],ff[MAXN],x[MAXN],xn[MAXN],xn 1[MAXN],f0,f1,f2,f3;for (i=0; i<N; i++) for (j=0; j<N; j++) if (j==i) ss[i][j]=1; else ss[i][j]=0;for (;;) {for (j=0; j<N; j++) xk[j]=x0[j];for (i=0; i<N; i++) {for (j=0; j<N; j++) sk[j]=ss[i][j];find_ab(0,1,&a,&b);search_gold(a,b,e2,&ax,&ay);for (j=0; j<N; j++) xk[j]=xkk[j];ff[i]=F(xk);}for (j=0; j<N; j++) xn[j] = xkk[j];for (j=0; j<N; j++) { sk[j]=xkk[j]-x0[j]; s1[j]=sk[j];}find_ab(0,1,&a,&b);search_gold(a,b,e2,&ax,&ay);for (j=0; j<N; j++) x[j]=xkk[j];d=0;for (j=0; j<N; j++) d+=(x[j]-x0[j])*(x[j]-x0[j]);d=sqrt(d);printf("k=%d;",k);for (j=0; j<N; j++) printf("x[%d]=%lf;",j+1,x0[j]);printf("d=%lf\n",d);if (d<=e1) {for (j=0; j<N; j++) x0[j]=x[j];break;}f0=F(x0);d=f0-ff[0]; m=0;for (j=1; j<N; j++) if (d<ff[j-1]-ff[j]) {m=j; d=ff[j-1]-ff[j];} for (j=0; j<N; j++) xn1[j]=2*xn[j]-x0[j];f1=F(x0); f2=F(xn); f3=F(xn1);if (0.5*(f1-2*f2+f3)>=d) {if (f2<f3) for (j=0; j<N; j++) x0[j]=xn[j];else for (j=0; j<N; j++) x0[j]=xn1[j];} else {for (i=m+1; i<N; i++) for (j=0; j<N; j++) ss[i-1][j]=ss[i][j];for (j=0; j<N; j++) ss[N-1][j]=s1[j];for (j=0; j<N; j++) x0[j]=x[j];}k++;}for (j=0; j<N; j++) x0[j]=xkk[j];return F(xkk);}main(){double fk,fkk,ck,d1,d2,e1,e2;double g[MAXN],h[MAXN],x1[MAXN],x0[MAXN];int i;N=2;nt=3;et=0;type=1;e1=0.0001;e2=0.0001;rk=0.001; ck=2;x0[0]=8; x0[1]=8;printf("=========\n");fk = nc_powell(x0,0.01,0.001);for(i=0; i<N; i++) x1[i]=x0[i];printf("rk=%lf\n",rk);for (i=0; i<N; i++) printf("x[%d]=%lf;",i+1,x0[i]);printf("%lf\n",fk);for(;;) {rk*=ck;fkk = nc_powell(x0,0.01,0.005);printf("rk=%lf\n",rk);for (i=0; i<N; i++) printf("x[%d]=%lf;",i+1,x0[i]); printf("%lf\n",fkk);d1=0; for(i=0; i<N; i++) d1+=(x1[i]-x0[i])*(x1[i]-x0[i]); d1=sqrt(d1);d2=fabs((fkk-fk)/fk);printf("d1=%lf d2=%lf\n",d1,d2);if(d1<=e1 || d2<=e2) break;for(i=0; i<N; i++) x1[i]=x0[i];fk=fkk;}fk=funt(x0,g,h);printf("**********************\n");for (i=0; i<N; i++) printf("x[%d]=%lf;",i+1,x0[i]);printf("F* = %lf\n",fk);for (i=0; i<nt; i++) printf("g[%d]=%lf;",i+1,g [i]); return 1;}。