鲍威尔方法求函数
鲍威尔法matlab实现
鲍威尔法% 鲍威尔法-计算第2环起始点和搜索方向syms xy1 xy2 fyx m x1 x2 X s1s2 Sf=60-10*x1-4*x2+x1^2+x2^2-x1*x2;disp ' 目标函数:f=60-10*x1-4*x2+x1^2+x2^2-x1*x2.'pretty(f);% 第1环沿e2方向搜索x01=0;x02=0;X0=[x01 x02]; % 第1环初始点x1=5.0000;x2=0;X=[x1 x2]; % 第1环沿e2方向初始点fy1=60-10*x1-4*x2+x1^2+x2^2-x1*x2; % 第1环e2方向初始点函数值e1=[1,0];e2=[0,1]; % 坐标轴单位方向s1=e2(1);s2=e2(2);S=[s1 s2];m=0;% 计算影射点及其函数值xs1=2*xy1-x01;xs2=2*xy2-x02;Xs=[xs1 xs2];fxs=60-10*xs1-4*xs2+xs1^2+xs2^2-xs1*xs2;fx0=60-10*x01-4*x02+x01^2+x02^2-x01*x02; % 第1环初始点函数值disp ' 第1环影射点坐标'disp (Xs)% 判断第1环搜索函数值下降最大方向df01=fx0-fy1; % 第1环沿e1方向初始点函数值-第1环沿e1方向终点函数值df02=fy1-fyx; % 第1环沿e2方向起始点函数值-第1环沿e2方向始点函数值if df01>df02dfm=df01;m=1;elsedfm=df02;m=2;end% 计算POWELL判别式f1=fx0;f2=fyx;f3=fxs;f123=[f1 f2 f3];disp ' 第1环初始点、终点、影射点函数值'disp (f123)fP1=(f1-2*f2+f3)*(f1-f2-dfm)^2;fP2=0.5*dfm*(f1-f3)^2;fP=[fP1 fP2];disp ' 两个POWELL判别式的值'fprintf(1,' POWELL判别式1的值 fP1= %10.4f \n',fP1)fprintf(1,' POWELL判别式1的值 fP2= %10.4f \n',fP2)fprintf(1,' 函数值下降最大方向 m= %2.0f \n',m)if f1>f3 & fP2>fP1sx1=xy1-X0(1); % 同时不满足两个判别式时56428 .产生新方向sx2=xy2-X0(2);sx=[sx1 sx2];if m==1S=[e2,sx]; % 新方向取代上环搜索中函数值下降最大方向m=1elseS=[e1,sx]; % 新方向取代上环搜索中函数值下降最大方向m=2 endelseS=[e1,e2]; % 满足某个判别式时56428 .保留上环搜索方向enddisp ' @@@@ 第2环搜索方向 @@@@'disp (S)% 第1环沿新方向搜索x1=xy1;x2=xy2;X=[x1 x2];s1=xy1-x01;s2=xy2-x02; S=[s1 s2];[xy1,xy2,fyx]=powell(m,x1,x2,X,s1,s2,S);disp '@@@@ 第2环搜索起始点 @@@@'if fyx<fxsx20=[xy1 xy2];disp ' (第1环极小点)'elsex20=[xs1 xs2];disp ' (第1环影射点)'enddisp (x20)% 鲍威尔法-第1环计算% 目标函数中常数项、步长的一次项和二次项系数fxs0=60-10*x1-4*x2+x1^2+x2^2-x1*x2;fxs1=-10*s1-4*s2+2*x1*s1+2*x2*s2-x1*s2-x2*s1;fxs2=s1^2+s2^2-s1*s2;fxs210=[fxs2 fxs1 fxs0];% 计算最优步长af=-fxs1/(2*fxs2);% 计算终点函数值(关于最优步长)fya=fxs0+fxs1*af+fxs2*af^2;% 计算终点坐标值和函数值(关于终点)xy1=x1+s1*af;xy2=x2+s2*af;Xopt=[xy1 xy2];fyx=60-10*xy1-4*xy2+xy1^2+xy2^2-xy1*xy2;if m==0disp ' ******** 计算第1环沿e2方向搜索 ********'elsedisp ' ******** 计算第1环沿新方向搜索 ********'enddisp ' 初始点'disp (X)disp ' 搜索方向'disp (S)disp ' 步长二次项一次项系数常数项'disp (fxs210)disp ' 最优步长'disp (af)if m==0disp ' 终点坐标'disp (Xopt)elsedisp ' 极小点坐标'disp (Xopt)disp ' 极小点函数值(关于最优步长)' disp (fya)disp ' 极小点函数值(关于极小点)' disp (fyx)end。
无约束优化方法程序
无约束优化方法---鲍威尔方法本实验用鲍威尔方法求函数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;}。
基于MATLAB的鲍威尔法求极值问题
基于MATLAB的鲍威尔法求极值问题：xxx 学号：xxx（理工大学机械与车辆学院车辆工程，100081）摘要：无约束优化方法主要有七种，按照求导与否把这些方法分为间接法和直接法。
牛顿法的成败与初始点选择有极大关系，其可靠性最差；坐标轮换法、单纯形法和最速下降法对于高维优化问题计算效率很低，有效性差；由于编制变尺度法程序复杂，其简便性不足。
综合考虑后，鲍威尔法、共轭梯度法具有较好的综合性能。
本文首先对鲍威尔法的原理进行阐述，根据其迭代过程给出流程图，并编写MATLAB程序。
最后用此MATLAB程序求解实际的极值问题，并对求解结果进行简要分析。
1.鲍威尔法的基本思想1.1其他优化方法对鲍威尔法形成的影响通过对鲍威尔法的学习，可以很明显看出来其迭代思想中汲取了其他几种优化方法的核心思想。
为了更全面、更深入的学习鲍威尔法，很有必要对其他有影响的优化思想进行学习和梳理。
由最基本的数学基础知识可知，梯度方向是函数增加最快的方向，负梯度方向是函数下降最快的方向，于是，利用这个下降最快方向产生了最速下降法。
每次迭代都沿着负梯度方向进行一维搜索，直到满足精度要求为止。
其特点是相邻两个搜索方向互相正交，所以很明显的一个现象就是刚开始搜索步长比较大，愈靠近极值点其步长愈小，收敛速度愈慢，特别当二维二次目标函数的等值线是较扁的椭圆时，迭代速度更慢。
这时，倘若目标函数是等值线长、短轴都平行于坐标轴的椭圆形，则通过坐标轮换法可以很高效的解决问题。
通过两次分别沿坐标轴进行一维搜索，便可达到极值点。
但对于目标函数的等值线椭圆的长、短轴倾斜于坐标轴时，坐标轮换法的搜索效率也显得极低。
抛开这两种特殊情况，对于一般形态的目标函数，如果在某些明显可以直达最优点的情况下（一般为靠近极值点区域），迭代过程完全可以不沿负梯度方向搜索，取而代之的是找到直达最优点的方向，一步到位。
但这样的直达方向应该如何去找呢？共轭梯度法由此产生。
其基本原理是：任意形式的目标函数在极值点附近的特性都近似一个二次函数，其等值线在极值点附近为近似的同心椭圆簇，而同心椭圆簇有一个特性便是任意两条平行线与椭圆簇切点的连线必通过椭圆的中心。
matlab实验鲍威尔法
实验报告实验名称：鲍威尔法院（系）：机电学院专业班级：机械制造及其自动化姓名：学号：2013年5 月13 日实验一：鲍威尔法实验日期：2013年5 月13 日一、实验目的了解MATLAB的基本运用了解MATLB在优化中的使用二、实验原理鲍威尔法也是一种共轭法，利用函数值来构造共轭方向，同时引入坐标轮换的概念，利用搜索前后两个点之间的连线形成新的共轭方向，替换旧的共轭方向。
三、实验内容鲍威尔法程序：x0=[12;10];xk=x0;ie=10^(-7);ae=1;%初始化搜索方向d=zeros(2,2);d(:,1)=[1;0];d(:,2)=[0;1];Inc=zeros(2,1);k=0;MLN=100;%迭代求解while (ae>ie&&k<MLN)syms x1syms x2xktemp=xk;fun1=fun(x1,x2);fun1=inline(fun1);f0=feval(fun1,xk(1),xk(2));F0=f0;if k>0F0=eval(F0);end%沿d1方向进行一维搜索syms asyms x1;syms x2;xk1=xk+a*d(:,1);x1=xk1(1);x2=xk1(2);fun1=fun(x1,x2);fxa=diff(fun1,'a');a=solve(fxa);xk1=inline(xk1);xk1=feval(xk1,a);xk1(1)=eval(xk1(1));xk1(2)=eval(xk1(2));syms x1;syms x2;fun1=fun(x1,x2);fun1=inline(fun1);f1=feval(fun1,xk1(1),xk1(2)); f1=eval(f1);Inc(1)=f0-f1;%沿d2方向进行搜索syms a;syms x1;syms x2;xk2=xk1+a*d(:,2);x1=xk2(1);x2=xk2(2);fun1=fun(x1,x2);fxa=diff(fun1,'a');a=solve(fxa);xk2=inline(xk2);xk2=feval(xk2,a);xk2(1)=eval(xk2(1));xk2(2)=eval(xk2(2));syms x1;syms x2;fun1=fun(x1,x2);fun1=inline(fun1);f2=feval(fun1,xk2(1),xk2(2));f2=eval(f2);F2=f2;Inc(2)=f1-f2;[Incm,row]=max(Inc);x3=2*xk2-xk;%计算反射点syms x1;syms x2;fun1=fun(x1,x2);fun1=inline(fun1);f3=feval(fun1,x3(1),x3(2));f3=eval(f3);F3=f3;temp1=(F0-2*F2+F3)*(F0-F2-Incm)^2; temp2=0.5*Incm*(F0-F3)^2;%判断是否更换搜索方向if (F3<F0&&temp1<temp2)syms a;syms x1;syms x2;d(:,row)=xk2-xk;xk=xk2+a*d(:,row);x1=xk(1);x2=xk(2);fun1=fun(x1,x2);fxa=diff(fun1,'a');a=solve(fxa);xk=inline(xk);xk=feval(xk,a);%不更换搜索方向else if F2<F3xk=xk2;elsexk=x3;endendxkerror=eval(xk2-xktemp); ae=norm(xkerror);k=k+1;endx=eval(xk)函数程序：function [f]=fun(x1,x2)f=2*x1^2+4*x1*x2+x2^2执行结果：x =四、实验小结通过本实验了解了了matlab的基本操作方法，了解鲍威尔法的原理与基本运用。
机械优化设计——鲍威尔法
机械优化设计——鲍威尔法机械优化设计是指在构筑机械产品或系统时，通过优化设计、算法分析等手段，使得机械产品或系统达到最佳的性能和效率。
鲍威尔法（Broyden-Fletcher-Goldfarb-Shanno algorithm，简称BFGS算法）是一种求解非线性最优化问题的拟牛顿法，广泛应用于机械优化设计中。
鲍威尔法是由四位科学家共同提出的，在非线性优化问题中具有很高的效率和可靠性。
它与传统的牛顿法相比，具有更好的收敛性能和更小的存储需求。
鲍威尔法是一种迭代的方法，通过不断更新设计变量的值，使得目标函数的值不断接近最优值。
鲍威尔法的核心思想是利用目标函数的梯度信息来指导方向和步长的选择。
在每一次迭代中，鲍威尔法通过构造一个正定的Hessian矩阵的近似来逼近目标函数的二次泰勒展开式，从而求取方向。
而步长的选择则通过线方法来确定，确保目标函数值在每次迭代中都能够下降。
在机械优化设计中，鲍威尔法可以应用于多种问题，如结构优化、参数优化等。
以结构优化为例，鲍威尔法可以通过不断调整结构的几何形状或内部结构的参数，来使得结构的重量最小或强度最大。
在进行参数优化时，鲍威尔法可以通过调整设计变量的取值，使得可以在给定约束条件下，最大化或最小化一些性能指标。
机械优化设计中的鲍威尔法需要满足一定的数学条件才能保证其收敛性和稳定性。
首先，目标函数必须是连续可微的。
其次，初始设计点的选取对于迭代过程的收敛性也有一定的影响。
一般情况下，选择一个合适的初始设计点可以加快鲍威尔法的收敛速度。
总结起来，鲍威尔法是机械优化设计中一种常用的求解非线性最优化问题的方法。
它可以通过迭代更新设计变量的值，使得目标函数达到最优值。
鲍威尔法不仅具有较好的收敛性能和计算效率，而且对于各种类型的机械优化设计问题都具有很好的适应性。
鲍威尔算法matlab程序-f
function f=fun(x)f=10*(x(1)+x(2)-5)^2+(x(1)-x(2))^2; function f=fx(x0,alpha,s)x1=x0+alpha*s;f=fun(x1);function f=fsearch(x0,s)%利用进退法确定高低高区间alpha1=0;h=0.1;alpha2=alpha1+h;f1=fx(x0,alpha1,s);f2=fx(x0,alpha2,s);if f1>f2alpha3=alpha2+h;f3=fx(x0,alpha3,s);while f2>f3alpha1=alpha2;alpha2=alpha3;alpha3=alpha3+h;f2=f3;f3=fx(x0,alpha3,s);endelseh=-h;v=alpha1;alpha1=alpha2; alpha2=v;v=f1;f1=f2;f2=v;alpha3=alpha2+h;f3=fx(x0,alpha3,s); while f2>f3alpha1=alpha2; alpha2=alpha3; alpha3=alpha3+h;f2=f3;f3=fx(x0,alpha3,s); endenda=min(alpha1,alpha3); b=max(alpha1,alpha3); %利用黄金分割点法求解alpha1=a+0.382*(b-a); alpha2=a+0.618*(b-a);f1=fx(x0,alpha1,s);f2=fx(x0,alpha2,s); while abs(a-b)>0.001 if f1>f2a=alpha1;alpha1=alpha2;f1=f2;alpha2=a+0.618*(b-a); f2=fx(x0,alpha2,s); elseb=alpha2;alpha2=alpha1;f2=f1;alpha1=a+0.382*(b-a); f1=fx(x0,alpha1,s); endendf=0.5*(a+b);clear%初始点x0=[0;0];%搜索方向e2=[0;1];G0=fun(x0);F0=G0;%第一次迭代%沿着e1alpha1=fsearch(x0,e1);x1=x0+alpha1*e1;F1=fun(x1);delta1=F0-F1;% 沿着方向e2；alpha2=fsearch(x1,e2);x2=x1+alpha2*e2;F2=fun(x2);G2=F2;delta2=F1-F2;deltam=max(delta1,delta2);%映射点x3=2*x2-x0;G3=fun(x3);if G3<G0 & (G0-2*G2+G3)*(G0-G2-deltam)^2<0.5*deltam*(G0-G3) ^2%方向替换e1=e2;e2=s;% 沿着方向s 进行搜索alpha3=fsearch(x2,s);x3=x2+alpha2*s;x0=x3;elseif F2>G3x0=x3;elsex0=x2;endEnd子文件JT,JH进退法程序代码�56555 .function [minx,maxx] = minJT(f,x0,h0,eps) format long;if nargin == 3eps = 1.0e-6;endk = 0;h = h0;while 1x4 = x1 + h;k = k+1;f4 = subs(f, findsym(f),x4); f1 = subs(f, findsym(f),x1); if f4 < f1x2 = x1;x1 = x4;f2 = f1;f1 = f4;h = 2*h;elseif k==1h = -h;x2 = x4;f2 = f4;elsex3 = x2;x2 = x1;break;endendendminx = min(x1,x3);maxx = x1+x3 - minx;format short;黄金分割法程序代码�56555 . function [x,minf] = minHJ(f,a,b,eps) format long;if nargin == 3eps = 1.0e-6;endl = a + 0.382*(b-a);u = a + 0.618*(b-a);k=1;tol = b-a;while tol>eps && k<100000fl = subs(f , findsym(f), l);fu = subs(f , findsym(f), u);a = l;l = u;u = a + 0.618*(b - a);elseb = u;u = l;l = a + 0.382*(b-a);endk = k+1;tol = abs(b - a);endif k == 100000disp('ÕÒ²»μ½215 ?î208 ?¡214 Oμ£¡); x = NaN;minf = NaN;return;endx = (a+b)/2;minf = subs(f, findsym(f),x);format short;。
004第四章Powell法
例4-2 74页
上述基本算法仅具有理论意义 :
3.改进的鲍威尔方法
在改进的算法中首先判断原向量组是否需要替换。 如果需要替换，还要进一步判断原向量组中 哪个向量最坏，然后再用新产生的向量替换这个 最坏的向量，以保证逐次生成共轭方向。
3.改进的鲍威尔方法
为此，要解决两个关键问题：
（1）dk＋1是否较好？是否应该进入新的方向组？ 即方向组是否进行更新？
0 1 2 2 x d d x 此轮基本方向组为 3 ， 3 ，起始点为 0 ＝ ，先 1 0 d 后沿 d 3 ， 3 方向，进行一维搜索，得
4 x ， 2
2 1
4 x 2
2 2
4 x ， 2
2 1
4 x 2
2 2
检验终止条件
x22 x02 0
3.96 x 1.9
1 2
f 2 f ( x1 ) 7.996 2
• （2）第2轮迭代计算
3.96 x 1.9
1 2
f 2 f ( x1 ) 7.996 2
构成新的方向
3.96 3.8 0.16 d x x 1.94 1.7 0.24
1.共轭方向
(d 0 )T Gd 1 0 (d 0 )T 2 f ( x)d 1 0
就是使d1直指极小点x* ， d1所必须满足的条件 。
两个向量 d 和d1称为G的共轭向量，或称
d 和d 对G是共轭方向。
0 1
0
共轭方向的性质
性质1 性质2 性质3 若非零向量系d0,d1,d2,…,dm-1是对G共轭， 在n维空间中互相共轭的非零向量的个数 从任意初始点出发，顺次沿n个G的共轭方 则这m个向量是线性无关的。
机械优化设计鲍威尔法
机械优化设计鲍威尔法
机械优化设计鲍威尔法（Powell method）是一种常用的非线性优化
算法，它以鲍威尔的名字命名，用于解决无约束非线性优化问题。
该方法
在各个领域都有广泛的应用，如工程优化设计、机器学习等。
下面将详细
介绍机械优化设计鲍威尔法的原理及应用。
鲍威尔法的具体步骤如下：
1.初始化参数：选择初始设计参数和方向。
2.寻找一维极小值点：沿着方向找到目标函数在该方向上的极小值点。
3.更新方向：通过比较前后两个极小值点的差异来更新方向。
4.迭代优化：重复步骤2和步骤3，直到达到指定的收敛条件。
鲍威尔法的优点是收敛速度较快、计算量较小，同时可以处理非线性
的优化问题。
然而，该方法也存在一些不足之处，如可能陷入局部最优解、对初值敏感等。
机械优化设计鲍威尔法在工程领域中有广泛的应用。
例如，在机械结
构设计中，可以利用鲍威尔法来优化结构参数，以满足特定的性能指标。
在汽车工业中，可以使用鲍威尔法来优化车辆的燃油效率和性能。
在航空
航天领域，可以利用该方法来优化飞行器的飞行性能。
此外，该方法还可
以用于机器学习中的参数优化，如调整神经网络的权重和偏置等。
总之，机械优化设计鲍威尔法是一种常用的非线性优化算法，通过迭
代逼近最优解。
虽然该方法有一些不足之处，但在实际应用中具有广泛的
适用性，尤其在工程优化设计和机器学习等领域。
通过使用该方法，可以
优化设计参数，改进性能指标，提高工程效率和产品质量。
鲍威尔法
根据这一原理构造的迭代算法称为鲍威尔基 本算法。
（二）鲍威尔法的缺陷
鲍威尔基本算法不可能对每一个都起作用，因为在迭代 过程中的n个搜索方向有时会变成线性相关的，而不能形 成共轭方向，导致随后的迭代搜索在降维（退化）的空间 中进行，可能求不到极小点，故而进行改进。
（三）鲍威尔修正算法
为了避免这种“退化”现象的发生，鲍威尔对这一算法 进行了修正。即在每一轮产生新的搜索方向 后，首先 判断原搜索方向组是否可以直接用下一轮迭代的方向组， 若可以，即用。否则，还要进一步判断原搜索方向组中哪 个方向上的函数值下降量最大，然后再用新搜索方向替换 这个下降量最大的搜索方向，以保证逐次生成共轭方向， 即每一轮迭代的搜索方向组线性无关。
并求出 （5）计算 判断
是否成立 若成立，则由 出发沿 方向进行一维搜索，求出目标 函数f(x)的极小点 ，并作为k+1轮的初始点 ，然后进行 k+1轮搜索，挤掉 ，同时把 放在方向组的最后 构成新一轮的方向组。
（6）若上述判断条件不成立，则k+1轮的初始点和方向组为 = 即此时k+1轮的n个搜索方向全部用第k轮的搜索方向。 （7）每轮迭代结束，都应检验收敛条件，若能满足：
鲍威尔法
鲍威尔（Powell）法又称方向加速度法， 它是利用共轭方向可以加快收敛速度的性质 形成的一种搜索方法。该方法不用对目标函 数求导，当目标函数的导数不连续时也能应 用，因此，鲍威尔法是一种十分有效的直接 搜索法。
一 共轭方向的概念与共轭向量的性质
（一）共轭方向 设A为n阶实对称正定矩阵，若有两个n维向量 和 能满足 A =0 则称向量 与 对矩阵A共轭，共轭向量的 方向称为共轭方向。
二 鲍威尔法
（一）鲍威尔法的基本原理和迭代过程
MATLAB鲍威尔算法
MATLAB鲍威尔算法鲍威尔算法（Broyden-Fletcher-Goldfarb-Shanno, BFGS）是用于非线性优化问题的一种数值优化算法。
它是一种拥有全局收敛性和快速收敛速度的准牛顿法。
BFGS算法的基本思想是通过近似二次函数来逼近原函数的局部结构，并利用此近似函数来求解极值。
它通过建立二次模型来估计目标函数的海森矩阵的逆（或近似逆），然后使用逆海森矩阵来更新方向。
算法的基本步骤如下：1.初始化参数：给定初始点x_0，设定精度要求ε，设置迭代次数k=0，以及初始H_0=I。
2.计算梯度：计算目标函数在当前点x_k的梯度g_k。
3.求解方向：计算方向d_k=-H_k*g_k，其中H_k表示当前的逆海森矩阵。
4.一维：在方向上进行一维线，求解最优步长α_k。
5.更新参数：更新参数x_{k+1}=x_k+α_k*d_k。
6.判断停止条件：如果，g_{k+1}，<ε，则停止迭代。
7. 更新逆海森矩阵：更新逆海森矩阵H_{k+1} = H_k + \frac{y_k* y_k^T}{y_k^T * s_k} - \frac{H_k * s_k * s_k^T * H_k}{s_k^T *H_k * s_k}，其中y_k = g_{k+1} - g_k，s_k = x_{k+1} - x_k。
8.如果迭代次数k超过预设迭代次数或者步长α_k小于预设步长，则停止迭代。
BFGS算法通过逆海森矩阵的更新来逼近目标函数的局部曲率，从而实现更快的收敛速度。
在每一次迭代中，BFGS算法更新逆海森矩阵，使其逐渐收敛到目标函数的真实海森矩阵的逆。
由于逆海森矩阵的计算复杂度为O(n^2)，其中n为变量的维度，BFGS算法在高维问题中的计算效率较低。
需要注意的是，BFGS算法只适用于无约束的非线性优化问题。
对于带有线性等式约束或者线性不等式约束的优化问题，需要使用相应的约束处理方法来进行求解。
总之，BFGS算法是一种拥有全局收敛性和快速收敛速度的准牛顿法。
机械优化设计鲍威尔法编程
}
}
// Set the icon for this dialog. The framework does this automatically
// when the application's main window is not a dialog
UpdateData(true);
int i,n;
double h1,h2,h3,m,flag,X00[2][1],d01[2][1],d02[2][1],d03[2][1]; //确定键的执行程序
double X01[2][1],X02[2][1],X03[2][1];
double F0,F1,F2,F3,e1,e2,em;
// the minimized window.
HCURSOR CMyDlg::OnQueryDragIcon()
{
return (HCURSOR) m_hIcon;
}
void CMyDlg::OnOK()
{
// TODO: Add extra validation here
//CDialog::OnOK();
F1=(X01[0][0]*X01[0][0]+2*X01[1][0]*X01[1][0]-4*X01[0][0]-2*X01[0][0]*X01[1][0]),f1=F1; //确定F的函数值
h2=(4*d02[0][0]+2*d02[0][0]*X01[1][0]+2*d02[1][0]*X01[0][0])/(2*d02[0][0]*d02[0][0]+ //确定搜索步长
十一、Powell算法（鲍威尔算法）原理以及实现
⼗⼀、Powell算法（鲍威尔算法）原理以及实现⼀、介绍 Powell算法是图像配准⾥⾯的常⽤的加速算法，可以加快搜索速度，⽽且对于低维函数的效果很好，所以本篇博客主要是为了介绍Powell算法的原理以及实现。
由于⽹上已经有了对于Powell算法的讲解，所以我只是把链接放出来（我觉得⾃⼰⽬前还没有这个讲解的能⼒），⼤家⾃⼰去了解。
放在这⾥主要也是为了节省⼤家搜索的时间。
（都是我⾟⾟苦苦搜出来的^-^）。
⼆、预备知识 了解⼀维搜索算法：进退法，消去法，黄⾦分割法 阅读以下博客：三、鲍威尔算法 具体原理阅读这⾥： 参考博客： 原理与例⼦（⼀个例⼦的计算过程）：四、matlab代码实现⼀个简单函数的求解 代码来源： 这个代码的程序与思路很是简洁，我觉得写得很好。
原⽂代码放在这⾥： ⽂件：MyPowell.m function MyPowell()syms x1 x2 x3 a;f=10*(x1+x2-5)^4+(x1-x2+x3)^2 +(x2+x3)^6;error=10^(-3);D=eye(3);x0=[000]';for k=1:1:10^6MaxLength=0;x00=x0;m=0;if k==1,s=D;endfor i=1:3x=x0+a*s(:,i);ff=subs(f,{x1,x2,x3},{x(1),x(2),x(3)});t=Divide(ff,a); %调⽤了进退法分割区间aa=OneDemensionslSearch(ff,a,t); %调⽤了0.618法进⾏⼀维搜索 xx=x0+aa*s(:,i);fx0=subs(f,{x1,x2,x3},{x0(1),x0(2),x0(3)});fxx=subs(f,{x1,x2,x3},{xx(1),xx(2),xx(3)});length=fx0-fxx;if length>MaxLength,MaxLength=length;m=m+1;endx0=xx;endss=x0-x00;ReflectX=2*x0-x00;f1=subs(f,{x1,x2,x3},{x00(1),x00(2),x00(3)});f2=subs(f,{x1,x2,x3},{x0(1),x0(2),x0(3)});f3=subs(f,{x1,x2,x3},{ReflectX(1),ReflectX(2),ReflectX(3)});if f3<f1&&(f1+f3-2*f2)*(f1-f2-MaxLength)^2<0.5*MaxLength*(f1-f3)^2x=x0+a*ss;ff=subs(f,{x1,x2,x3},{x(1),x(2),x(3)});t=Divide(ff,a);aa=OneDemensionslSearch(ff,a,t);x0=x0+aa*ss;for j=m:(3-1),s(:,j)=s(:,j+1);ends(:,3)=ss;elseif f2>f3, x0=ReflectX;endendif norm(x00-x0)<error,break;endk;x0;endopx=x0;val=subs(f,{x1,x2,x3},{opx(1),opx(2),opx(3)});disp('最优点：');opx'disp('最优化值：');valdisp('迭代次数：');k ⽂件 Divide.m ：%对任意⼀个⼀维函数函数进⾏区间分割，使其出现“⾼—低—⾼”的型式function output=Divide(f,x,m,n)if nargin<4,n=1e-6;endif nargin<3,m=0;endstep=n;t0=m;ft0=subs(f,{x},{t0});t1=t0+step;ft1=subs(f,{x},{t1});if ft0>=ft1t2=t1+step;ft2=subs(f,{x},{t2});while ft1>ft2t0=t1;%ft0=ft1;t1=t2;ft1=ft2;step=2*step;t2=t1+step;ft2=subs(f,{x},{t2});endelsestep=-step;t=t0;t0=t1;t1=t;ft=ft0;%ft0=ft1;ft1=ft;t2=t1+step;ft2=subs(f,{x},{t2});while ft1>ft2t0=t1;%ft0=ft1;t1=t2;ft1=ft2;step=2*step;t2=t1+step;ft2=subs(f,{x},{t2});endendoutput=[t0,t2];View Code ⽂件：OneDemensionslSearch.mfunction output=OneDemensionslSearch(f,x,s,r)if nargin<4,r=1e-6;enda=s(1);b=s(2);a1=a+0.382*(b-a);fa1=subs(f,{x},{a1});a2=a+0.618*(b-a);fa2=subs(f,{x},{a2});while abs((b-a)/b)>r && abs((fa2-fa1)/fa2)>rif fa1<fa2b=a2;a2=a1;fa2=fa1;a1=a+0.382*(b-a);fa1=subs(f,{x},{a1});elsea=a1;a1=a2;fa1=fa2;a2=a+0.618*(b-a);fa2=subs(f,{x},{a2});endendop=(a+b)/2;%fop=subs(f,{x},{op});output=op;View Code 全部放到同⼀个⼯程⽬录⾥⾯，设置为当前⽬录，然后输⼊Powell即可运⾏得到结果。
机械优化设计方法——鲍威尔法,坐标轮换法
鲍威尔法
机械优化设计方法——鲍威尔法
鲍威尔(Powell)法是直接利用函数值 来构造共扼方向的一种共扼方向法。基本 思想是在不用求导数的前提下，在迭代中 逐次构造Hessian矩阵G的共扼方向。
原始共轭法的缺点
对于n维无约束最优化问题，采用原始共轭 方向法在产生n个共轭方向时，有可能是线 性相关或接近线性相关的，遇这种情况， 会导致在降维空间寻优使迭代计算不能收 敛到真正最优点而失败 。
(4)判断是否满足收敛准则。
即
X k 1 X k
。若满足，则
X k1 0
为
极小点，否则应置 k k 1 ，返回，
继续进行下一轮迭代。
计算 f (X ) x12 25x22 的极小值
解 ： 选取初始点为 X 0 2,2T,初试搜索方向
为 d10 1,0T，d20 0,1 T。初始点处 的函数
值
F1
f
(
X
0 0
)
104
第一轮迭代：沿 d10 方向进行一维搜
索，得
X10
X
0
1d10
2
1
1 0
2
1
2
将
X
0 1
代入原方程得：
f1 f (X10) (2 1)2 2522 4 41 12 100
最佳步长
：
df ( d
)
4
21
0
1 2
X10
2 2
210
0 2
X
0 1
鲍威尔法计算步骤与程序原理
(1)任选一初始点 X 0 ,
X
0 0
X
0
，选取初
试方向组，它由n个线性无关的向量
matlab鲍威尔法求二元二次函数的极小值
matlab鲍威尔法求二元二次函数的极小值鲍威尔法（Powell's method）是一种用于求解无约束优化问题的迭代算法。
在MATLAB中，你可以使用内建函数，比如fminunc 或fminsearch，或者手动实现鲍威尔法来求解二元二次函数的极小值。
不过，MATLAB并没有直接提供鲍威尔法的内建函数，因此你需要自己实现它。
下面是一个简化的鲍威尔法实现，用于求解二元二次函数的极小值。
假设我们的二元二次函数是f(x, y) = ax^2 + bxy + cy^2 + dx + ey + f。
matlabfunction [xmin, fmin] = powell_method(func, x0, tol, max_iter) % func: 目标函数句柄% x0: 初始点% tol: 收敛容忍度% max_iter: 最大迭代次数n = length(x0); % 变量的维数x = x0;B = eye(n); % 初始基矩阵为单位矩阵for k = 1:max_iter% 在当前基方向上进行一维搜索for i = 1:n% 定义搜索方向d = B(:, i);% 一维线搜索确定步长alpha = linesearch(func, x, d);% 更新当前点x_new = x + alpha * d;% 检查收敛性if norm(x_new - x) < tolreturnend% 更新xx = x_new;% 更新基矩阵B(:, n) = x_new - x;B = B(:, 1:n-1);end% 使用QR分解更新基矩阵[Q, R] = qr(B);B = Q(:, 1:n);endxmin = x;fmin = func(x);endfunction alpha = linesearch(func, x, d)% 简单的线搜索实现（这里假设函数是凸的）alpha = 0.1; % 初始步长c1 = 1e-4; % 足够小的正数while func(x + alpha * d) > func(x) + c1 * alpha * d' * grad(func, x)alpha = alpha / 2;endendfunction g = grad(func, x)% 计算梯度（这里需要func的梯度信息）% 对于二次函数ax^2 + bxy + cy^2 + dx + ey + f% 梯度是[2ax + by + d, bx + 2cy + e]'% 这里只是一个示例，你需要根据实际的func来计算梯度% 假设a = 1;b = 2;c = 1;d = -4;e = -6;g = [2*a*x(1) + b*x(2) + d; b*x(1) + 2*c*x(2) + e];end% 示例：求解函数f(x, y) = x^2 + 2xy + y^2 - 4x - 6y 的极小值func = @(x) x(1)^2 + 2*x(1)*x(2) + x(2)^2 - 4*x(1) - 6*x(2);x0 = [0; 0]; % 初始点tol = 1e-6; % 容忍度max_iter = 1000; % 最大迭代次数% 调用鲍威尔法[xmin, fmin] = powell_method(func, x0, tol, max_iter);% 显示结果disp(['极小值点:', mat2str(xmin)]);disp(['函数极小值:', num2str(fmin)]);请注意，上面的代码片段中有几个地方需要特别注意：grad 函数需要根据你的目标函数来计算梯度。
改进鲍威尔法
x=[1;1];%初始点
ff(1)=f(x(1),x(2));
e=0.001;%精度为0.001
d=[1;0;0;1];
while(1)
x00=[x(1);x(2)];
fori=1:n
[a(i),b(i)]=section(x(2*i-1),x(2*i),d(2*i-1),d(2*i));
alpha(i)=goldencut(x(2*i-1),x(2*i),d(2*i-1),d(2*i),a(i),b(i));
fori=1:n
ifdelta==Delta(i)
m=i;
break;
end
end
d(2*n+1)=x(2*n+1)-x(1);
d(2*n+2)=x(2*n+2)-x(2);
x(2*n+3)=2*x(2*n+1)-x(1);
x(2*n+4)=2*x(2*n+2)-x(2);
ff(n+2)=f(x(2*n+3),x(2*n+4));
f0=ff(1);f2=ff(n+1);f3=ff(n+2);
k=k+1;
R(k,:)=[k,x',d',ff];
If f3<f0 &(f0-2*f2+f3)*(f0-f2-delta)^2<0.5*delta*(f0-f3)^2
[a(n+1),b(n+1)]=section(x(2*n+1),x(2*n+2),d(2*n+1),d(2*n+2));
alpha2=a+r*(b-a);
鲍威尔法概述及算例求解
一 共轭方向的概念与共轭向量的性质
（一）共轭方向 设A为n阶实对称正定矩阵，若有两个n维向量 和 能满足 A =0 则称向量 与 对矩阵A共轭，共轭向量的 方向称为共轭方向。
（二）共轭向量的性质 设A为n×n阶实对称正定矩阵， （i=1,2,…n） 是关于A的n个互相共轭的非零向量，对于正 定二次函数f(x)的极小化寻优问题，从任何初 始点出发，依次沿 方向经n次一维搜索即可 收敛到极小点 = 沿n元二次正定函数的n个共轭方向进行n次 一维搜索就可达到目标函数的极小点。
则可输出最优解，迭代结束。否则进入下一轮迭代，即转 步骤（2）
（5）计算举例 例：用鲍威尔法求目标函数 解（极小值）。已知，初始点
的最优 ，收敛精度
满足
进行第二次循环
ห้องสมุดไป่ตู้
总结 ： （1）共轭方向及共轭向量 （2）鲍威尔修正算法，判别条件及计算步骤。
n二鲍威尔法一鲍威尔法的基本原理和迭代过程1采用坐标轮换法顺次沿n个坐标方向进行一维搜索然后以初始点和终点构成一个新的方向方向为搜索方向再作一维搜索得到极小点2取始点并去掉原搜索方向组中的第一个方向而将第一轮构成的新搜索方向方向而将第轮构成的新搜索方向末一个方向以此组成第二轮迭代的n个方向
鲍威尔法
鲍威尔（Powell）法又称方向加速度法， 它是利用共轭方向可以加快收敛速度的性质 形成的一种搜索方法。该方法不用对目标函 数求导，当目标函数的导数不连续时也能应 用，因此，鲍威尔法是一种十分有效的直接 搜索法。
根据这一原理构造的迭代算法称为鲍威尔基 本算法。
（二）鲍威尔法的缺陷
鲍威尔基本算法不可能对每一个都起作用，因为在迭代 过程中的n个搜索方向有时会变成线性相关的，而不能形 成共轭方向，导致随后的迭代搜索在降维（退化）的空间 中进行，可能求不到极小点，故而进行改进。
MATLAB关于Powell优化设计程序
f0 = X0_1(1)^2+2*X0_1(2)^2-4*X0_1(1)-2*X0_1(1)*X0_1(2); % 初始点的函数值.
Dt1_1 = f0-f1; % 计算本轮相邻两点函数值的下降量. Dt2_1 = f1-f2; Dtm_1 = max(Dt1_1,Dt2_1); % 进行收敛判断(是否用得到的第一个共轭方向替换上一轮中的第一个一维搜索方向). if (F3<F1&&(F1+F3-2*F2)*(F1-F2-Dtm_1)^2<0.5*Dtm_1*(F1-F3)^2) S1_2 = S2_1; S2_2 = S_1; else S1_2 = S2_1; S2_2 = S1_1; end syms a3 % 以下语句是求出沿S_1方向进行一维搜索的最佳步长因子以及第二轮迭代的初始点X0_2. X_1 = X2_1+a3*S_1; f3 = X_1(1)^2+2*X_1(2)^2-4*X_1(1)-2*X_1(1)*X_1(2); ff3 = diff(f3); a3 = solve(ff3,0); % 求得沿S_1方向进行一维搜索的最优步长 a3. X_1 = X2_1+a3*S_1; f3 = eval(X_1(1)^2+2*X_1(2)^2-4*X_1(1)-2*X_1(1)*X_1(2)); % 得到第二轮迭代的初始点X_1处的函数值. X0_2 =eval(X_1); F_1 = f3; % 进行迭代终止检验 d1 = sqrt((X0_2(1)-X0_1(1))^2+(X0_2(2)-X0_1(2))^2); if (d1>E) % 得到d1 = 2.886173937932362 fprintf('第一轮迭代完成过后的精度检验值为：d1 = %4f\n',d1) % 进行迭代终止检验是否继续进行下一轮迭代 % 第二轮迭代(K=2!) % 沿S2_1方向进行第二轮迭代第一次一维搜索 syms a4 % 以下语句是求出沿S1_2方向进行一维搜索的最佳步长因子 X1_2 = X0_2+a4*S1_2; f4 = X1_2(1)^2+2*X1_2(2)^2-4*X1_2(1)-2*X1_2(1)*X1_2(2); ff4 = diff(f4); a4 = solve(ff4,0);% 得到第二轮迭代第一次一维搜索的最优步长因子a4. fprintf('第二轮迭代第一次一维搜索的最优步长因子为: a4 = %4f\n',eval(a4)) X1_2 = X0_2+a4*S1_2; f4 = eval(X1_2(1)^2+2*X1_2(2)^2-4*X1_2(1)-2*X1_2(1)*X1_2(2)); % 得到第二轮迭代点X1_2处的函数值f4. % 沿S2_2方向进行第二轮迭代第二次一维搜索 syms a5 X2_2 = X1_2 + a5*S2_2; f5 = X2_2(1)^2+2*X2_2(2)^2-4*X2_2(1)-2*X2_2(1)*X2_2(2); ff5 = diff(f5); % 得到第二轮的迭代初始点X0_2.
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
本实验用鲍威尔方法求函数f(x)=(x1-5)2+(x2-6)2 的最优解#include <math.h>
#include <stdio.h>
#include <stdlib.h>
const MAXN = 10;
double xkk[MAXN],xk[MAXN],sk[MAXN];
int N;
double F(double *x)
{
return 4*pow(x[0]-5,2.0)+pow(x[1]-6,2.0);
}
double f(double x)
{
for (int i=0; i<N; i++) xkk[i]=xk[i]+x*sk[i];
return F(xkk);
}
/*
无约束坐标轮换法
x0--初始点
e1--一维搜索精度
e2--求解精度
*/
double nc_trans(double *x0,double e1,double e2)
{
int i,j,k=1;
double a,b,ax,ay,d;
for (;;) {
for (j=0; j<N; j++) xk[j]=x0[j];
for (i=0; i<N; i++) {
for (j=0; j<N; j++) if (j==i) sk[j]=1;
else sk[j]=0;
find_ab(0,1,&a,&b);
search_gold(a,b,e2,&ax,&ay);
for (j=0; j<N; j++) xk[j]=xkk[j];
}
d=0;
for (j=0; j<N; j++)
d+=(x0[j]-xkk[j])*(x0[j]-xkk[j]);
d=sqrt(d);
printf("k=%d;",k);
for (j=0; j<N; j++)
printf(",x[%d]=%lf;",j+1,xkk[j]);
printf("d=%lf\n",d);
if (d<=e1) break;
for (j=0; j<N; j++) x0[j]=xkk[j];
k++;
}
for (j=0; j<N; j++) x0[j]=xkk[j];
return F(xkk);
}
/*
鲍威尔法
x0--初始点
e1--一维搜索精度
e2--求解精度
*/
double nc_powell(double *x0,double e1,double e2) {
int i,j,k=1,m;
double a,b,ax,ay,d;
double ss[MAXN][MAXN],s1[MAXN],
ff[MAXN],x[MAXN],xn[MAXN],
xn1[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);
}。