机械优化设计powell法程序
机械优化设计powell法程序
#include"stdio.h"#include"stdlib.h"#include"math.h"double objf(double x[]){double ff;ff=x[0]*x[0]+2*x[1]*x[1]-2*x[0]*x[1]-4*x[0];return(ff);}double gold(double a[],double b[],double eps,int n,double xx[]) {int i;double f1,f2,*x[2],ff,q,w;for(i=0;i<2;i++)x[i]=(double *)malloc(n*sizeof(double));for(i=0;i<n;i++){*(x[0]+i)=a[i]+0.618*(b[i]-a[i]);*(x[1]+i)=a[i]+0.382*(b[i]-a[i]);}f1=objf(x[0]);f2=objf(x[1]);do{if(f1>f2){for(i=0;i<n;i++){b[i]=*(x[0]+i);*(x[0]+i)=*(x[1]+i);}f1=f2;for(i=0;i<n;i++)*(x[1]+i)=a[i]+0.382*(b[i]-a[i]);f2=objf(x[1]);}else{ for(i=0;i<n;i++){a[i]=*(x[1]+i);*(x[1]+i)=*(x[0]+i);}f2=f1;for(i=0;i<n;i++)*(x[0]+i)=a[i]+0.618*(b[i]-a[i]);f1=objf(x[0]);}q=0;for(i=0;i<n;i++)q=q+(b[i]-a[i])*(b[i]-a[i]);}while(w>eps);for(i=0;i<n;i++)xx[i]=0.5*(a[i]+b[i]);ff=objf(xx);for(i=0;i<2;i++)free(x[i]);return(ff);}void jtf(double x0[],double h0,double s[],int n,double a[],double b[]) {int i;double *x[3],h,f1,f2,f3;for(i=0;i<3;i++)x[i]=(double *)malloc(n*sizeof(double));h=h0;for(i=0;i<n;i++)*(x[0]+i)=x0[i];f1=objf(x[0]);for(i=0;i<n;i++)*(x[1]+i)=*(x[0]+i)+h*s[i];f2=objf(x[1]);if(f2>=f1){h=-h0;for(i=0;i<n;i++)*(x[2]+i)=*(x[0]+i);f3=f1;for(i=0;i<n;i++){*(x[0]+i)=*(x[1]+i);*(x[1]+i)=*(x[2]+i);}f1=f2;f2=f3;}for(;;){h=2*h;for(i=0;i<n;i++)*(x[2]+i)=*(x[1]+i)+h*s[i];f3=objf(x[2]);if(f2<f3) break;else{ for(i=0;i<n;i++){*(x[0]+i)=*(x[1]+i);*(x[1]+i)=*(x[2]+i);f1=f2;f2=f3;}}if(h<0)for(i=0;i<n;i++){a[i]=*(x[2]+i);b[i]=*(x[0]+i);}elsefor(i=0;i<n;i++){a[i]=*(x[0]+i);b[i]=*(x[2]+i);}for(i=0;i<3;i++)free(x[i]);}double oneoptim(double x0[],double s[],double h0,double epsg,int n,double x[]) {double *a,*b,ff;a=(double *)malloc(n*sizeof(double));b=(double *)malloc(n*sizeof(double));jtf(x0,h0,s,n,a,b);ff=gold(a,b,epsg,n,x);free(a);free(b);return (ff);}double powell(double p[],double h0,double eps,double epsg,int n,double x[]) {int i,j,m;double *xx[4],*ss,*s;double f,f0,f1,f2,f3,fx,dlt,df,sdx,q,d;ss=(double *)malloc(n*(n+1)*sizeof(double));s=(double *)malloc(n*sizeof(double));for(i=0;i<n;i++){for(j=0;j<=n;j++)*(ss+i*(n+1)+j)=0;*(ss+i*(n+1)+i)=1;}for(i=0;i<4;i++)xx[i]=(double *)malloc(n*sizeof(double));for(i=0;i<n;i++)*(xx[0]+i)=p[i];{for(i=0;i<n;i++){*(xx[1]+i)=*(xx[0]+i);x[i]=*(xx[1]+i);}f0=f1=objf(x);dlt=-1;for(j=0;j<n;j++){for(i=0;i<n;i++){*(xx[0]+i)=x[i];*(s+i)=*(ss+i*(n+1)+j);}f=oneoptim(xx[0],s,h0,epsg,n,x);df=f0-f;if(df>dlt){dlt=df;m=j;}}sdx=0;for(i=0;i<n;i++)sdx=sdx+fabs(x[i]-(*(xx[1]+i)));if(sdx<eps){free(ss);free(s);for(i=0;i<4;i++)free(xx[i]);return(f);}for(i=0;i<n;i++)*(xx[2]+i)=x[i];f2=f;for(i=0;i<n;i++){*(xx[3]+i)=2*(*(xx[2]+i)-(*(xx[1]+i))); x[i]=*(xx[3]+i);}fx=objf(x);f3=fx;q=(f1-2*f2+f3)*(f1-f2-dlt)*(f1-f2-dlt); d=0.5*dlt*(f1-f3)*(f1-f3);if((f3<f1)||(q<d)){if(f2<=f3)for(i=0;i<n;i++)*(xx[0]+i)=*(xx[2]+i);for(i=0;i<n;i++)*(xx[0]+i)=*(xx[3]+i);}else{for(i=0;i<n;i++){*(ss+(i+1)*(n+1))=x[i]-(*(xx[1]+i));*(s+i)=*(ss+(i+1)*(n+1));}f=oneoptim(xx[0],s,h0,epsg,n,x);for(i=0;i<n;i++)*(xx[0]+i)=x[i];for(j=m+1;j<=n;j++)for(i=0;i<n;i++)*(ss+i*(n+1)+j-1)=*(ss+i*(n+1)+j);}}}void main(){double p[]={1,2};double ff,x[2];ff=powell(p,0.3,0.001,0.0001,2,x);printf("x[0]=%f,x[1]=%f,ff=%f\n",x[0],x[1],ff); getchar();}。
机械优化设计鲍威尔程序
ff=fny(xx); printf("\n\nThe Optimal Design Result Is:\n"); for(i=0;i<n;i++) printf("\n\tx[%d]*=%f",i+1,xx[i]); printf("\n\tf*=%f",ff); getch(); }
d=(pow(a[0],2)-pow(a[2],2))*(a[0]-a[1])-(pow(a[0],2)-pow(a[1],2))*(a[0]-a[2]); if(fabs(d)==0) break; c=((f[0]-f[2])*(a[0]-a[1])-(f[0]-f[1])*(a[0]-a[2]))/d; if(fabs(c)==0) break; b=((f[0]-f[1])-c*(pow(a[0],2)-pow(a[1],2)))/(a[0]-a[1]); aa=-b/(2*c); *ft=func(xk,aa,s); if(fabs(aa-a[1])<=ad) { if(*ft>f[1]) aa=a[1]; break; } if(aa>a[1]) { if(*ft>f[1]) { a[2]=aa; f[2]=*ft; } else if(*ft<f[1]) {a[0]=a[1]; a[1]=aa; f[0]=f[1]; f[1]=*ft; } else if(*ft==f[1]) { a[2]=aa; a[0]=a[1]; f[2]=*ft; f[0]=f[1]; a[1]=(a[0]+a[2])/2; f[1]=func(xk, a[1],s); } } else { if(*ft>f[1]) { a[0]=aa; f[0]=*ft; } else if (*ft<f[1]) { a[2]=a[1]; a[1]=aa; f[2]=f[1]; f[1]=*ft; } else if (*ft==f[1]) { a[0]=aa; a[2]=a[1]; f[0]=*ft; f[2]=f[1]; a[1]=(a[0]+a[2])/2; f[1]=func(xk,a[1],s); } }}
机械优化设计鲍威尔法编程
机械优化设计鲍威尔法编程鲍威尔法(Powell's method)是一种常用于机械优化设计的迭代算法,它基于步长的方向进行,进而找到局部或全局最优解。
该算法主要用于解决无约束优化问题,即不涉及约束条件的优化设计。
下面将详细介绍鲍威尔法的编程实现。
鲍威尔法的基本思路是在迭代过程中通过多次步长方向,找到全局最优解。
具体步骤如下:1.初始化:设置初始点x0和迭代次数k=0。
2.计算方向:选择一个初始的方向d0和步长α,并将d0归一化为单位向量。
3. 求解新的迭代点:通过计算当前点xk加上步长α乘以方向dk,得到新的迭代点xk+14. 更新方向:计算新的方向dk+15. 判断是否达到终止条件:如果达到了终止条件,则输出当前点xk+1为最优解;否则,令k=k+1,返回第3步继续进行迭代。
下面给出一个使用Python编程实现鲍威尔法的示例代码:```pythonimport numpy as npdef powell_method(f, x0, alpha, eps, max_iter):#初始化x=x0d = np.eye(len(x0))k=0while k < max_iter:#计算方向和步长g=f(x)d_norm = np.linalg.norm(d, axis=0) d = d / d_normalpha = alpha / d_norm#求解新的迭代点x_new = x + alpha * d#更新方向g_new = f(x_new)delta = g_new - gd = np.roll(d, -1, axis=0)d[-1] = (x_new - x) / alpha#判断终止条件if np.linalg.norm(delta) < eps: return x_new#更新迭代点x = x_newk+=1return x#示例函数,目标是求解f(x)=(x[0]-1)^2+(x[1]-2)^2 def f(x):return (x[0] - 1) ** 2 + (x[1] - 2) ** 2#设置初始点、步长、终止条件和最大迭代次数x0 = np.array([0.0, 0.0])alpha = 0.1eps = 1e-6max_iter = 100#调用鲍威尔法进行优化设计x_opt = powell_method(f, x0, alpha, eps, max_iter) #输出最优解print("Optimal solution: ", x_opt)print("Optimal value: ", f(x_opt))```在上述代码中,目标函数f(x)为示例函数,可以根据具体的优化设计问题进行修改。
优化设计-(Powell)法+鲍威尔修正算法
1)任选初始点X(0)=X0(1) ,给定迭代收敛精度 1 , 2 。 取初始基本方向组为单位坐标向量系,即Si(1)=ei ( i = 1, 2, … , n ), 并置迭代轮次k=1。
2)从X0(k)出发,依次沿Si(k) ( i = 1, 2, … , n ) 作一维 搜索,得n个极小点Xi(k) ( i = 1, 2, … , n ),构造新的搜 索方向 S(k) = Xn(k) -X0(k) ,并沿此方向进行一维搜索得 极小点Xn+1(k) 。 3)判断迭代终止条件 || Xn+1(k) – X0(k) ||≤ 1 ? 或 | f (Xn+1(k) ) – f (X0(k) ) |≤ 2 | f (Xn+1(k) ) | ? 若满足,则终止迭代并输出最优解: X * = Xn+1(k) 和 f * = f (X * )
依此进行下去,直到获得满足迭代收敛精度要求的 近似极小点为止。 根据这一原理构造的迭代算法称为鲍威尔基本 算法。 x2
S(1) S2(1) X3 (1) X2 (1) X2
* X (2)
S(2)
X1 (2) S1(1)
()
X0
(1)
X1
x1
2.3.3鲍威尔基本算法的缺点 鲍威尔基本算法仅具有理论意义,不要说对于 一般的函数,就是对于二次函数,它也可能失效。 因为在迭代过程中的 n个搜索方向有时会变成线性相
X0(k+1) = Xn+1(k) Si(k+1) ( i = 1, 2, … , n ) 即 { S1(k), …, Sm-1(k), S(k), Sm+1(k), …, Sn(k)} 并置 k k+1 ,返回第2)步。
机械优化设计.第十章
X (k) H (k 1) g (k) _ 拟牛顿条件(" DFP " 条件)
三、 ”DFP”变尺度法
17
1、 ”DFP”变尺度矩阵的递推公式
H (k 1) H (k ) △H (k )
(k 0,1,2,...)
X (k) H (k 1) g (k) _ 拟牛顿条件(" DFP " 条件)
5 无约束优化方法 5-3 powell法 一、共轭方向 二、Powell 基本算法 三、改进后的Powell算法
1
5 无约束优化方法
5-4 梯度法、牛顿法
5-5 DFP变尺度法简介
2
1、熟悉无约束优化方法梯度法与牛顿法的 基本思路及分类
2、比较间接法与直接法的特点及实用范围
3、了解变尺度法的基本思想,能够用其解 决无约束优化问题
二次函数(X)近似 代替原目标函数f ( X )
X
()
为f ( X )的下一个迭代新 点
X (k 1)
如此反复迭代,逐步逼近f ( X )的最优点X
8
2、迭代方法
二阶Taylor展开式:
f (X)
(X)
f (X (k) )
f (X (k) ) T X
1 2
X
T2
f
(X
(k ) )X
X [ x1,x2, xn ]T
X (k ) (H (k ) H (k ) )g (k ) H (k ) g (k ) X (k ) H (k ) g (k ) (1)
H (k) X (k) [q (k) ]T H (k) g (k) [w (k) ]T (2)
powel法
powel法Powell法是一种用于无约束优化问题的迭代算法,通过不断地寻找搜索方向和步长来逐步逼近最优解。
该算法在实际应用中具有广泛的适用性和高效性,因此被广泛应用于工程、经济、物理等领域。
其基本思想是将多个搜索方向组合起来构成一个新的搜索方向,从而使得每次迭代可以更加精确地逼近最优解。
具体来说,Powell法采用共轭方向法或方向加速法,通过沿连接相邻两轮搜索末端的向量S方向搜索,以共轭方向打破振荡,加速收敛。
在每轮迭代中,先沿初始方向组Si(1) (i=1,2,…,n)的n个方向和共轭方向S(1),搜索n+1次得极值点xn+1(1);然后沿方向组Si(2) ( i=1,2,…,n;i≠m )的n-1个方向和共轭方向S(1),构筑共轭方向S(2)搜索n+1次得极值点xn+1(2)。
Powell法的收敛速度一般较快,适用于求解大规模的无约束优化问题。
Powell法的基本步骤如下:1. 初始化:选择初始点x0,初始方向组Si(1) (i=1,2,…,n),以及初始步长λ1。
2. 迭代:对于k=1,2,…,进行以下步骤:a. 沿初始方向组Si(k)的n个方向和共轭方向Sk搜索n+1次,得到极值点xk+1(1)。
b. 计算搜索方向组Si(k+1) (i=1,2,…,n;i≠m),其中m为使Sk为最小值的下标。
c. 沿共轭方向Sk+1搜索n+1次,得到极值点xk+1(2)。
3. 判断停止条件:检查是否满足停止条件,例如达到预设的最大迭代次数或相邻两次迭代之间的差距小于预设的阈值。
如果满足停止条件,则返回当前点作为最优解;否则,继续迭代。
4. 更新:根据当前极值点和搜索方向更新λk+1和Sk+1,然后回到步骤2。
Powell法的优点在于其收敛速度较快,适用于求解大规模的无约束优化问题。
然而,该算法对初始点选择较为敏感,如果初始点选择不当,可能会导致算法无法收敛到全局最优解。
此外,Powell法对于某些特定的问题可能需要调整方向组的选择和搜索步长的确定方式。
powell法matlab程序
powell法matlab程序Powell法是一种用于求解无约束最优化问题的迭代优化算法。
它通过逐步旋转坐标轴的方式来寻找函数的最小值点。
在本文中,我们将详细介绍Powell法的原理和应用,并提供MATLAB程序实现。
Powell法的基本思想是通过旋转坐标轴的方式,将多维优化问题转化为一维优化问题。
这种方法通过变换坐标轴,将迭代过程中的每次更新只涉及到一个变量,从而降低了计算的复杂性。
Powell法在某些问题上比前一种单纯形装囊法更加高效。
下面是我们使用MATLAB实现Powell法的步骤:步骤1:定义目标函数首先,我们需要定义目标函数。
目标函数可以是任何连续可导的函数。
在MATLAB中,我们可以通过函数句柄(即指向目标函数的指针)来表示目标函数。
例如,我们可以定义一个简单的目标函数如下:matlabfunction f = myFunction(x)目标函数为x^2 + 2*x + 1f = x.^2 + 2.*x + 1;end步骤2:初始化参数接下来,我们需要初始化Powell法的参数。
这些参数包括初始点的位置和搜索方向。
我们可以选择任意合适的初始点,以及初始搜索方向。
例如,我们可以将初始点的位置设置为(1, 1)。
matlab初始点的位置x0 = [1, 1];初始搜索方向d = [1, 0];步骤3:定义搜索函数我们还需要定义一个搜索函数,用于根据当前位置和搜索方向来计算最佳的步长。
在Powell法中,可以使用一维搜索方法来寻找步长。
这里,我们可以使用黄金分割法(golden section method)来实现。
matlabfunction [alpha, val] = lineSearch(x, d)黄金分割法的参数rho = (sqrt(5) - 1) / 2;epsilon = 1e-6;定义目标函数f = (t) myFunction(x + t.*d);在初始区间上进行黄金分割法搜索a = 0;b = 1;h = b - a;c = a + rho.*h;d = b - rho.*h;while (h > epsilon)if (f(c) < f(d))b = d;elsea = c;endh = b - a;c = a + rho.*h;d = b - rho.*h;endalpha = (a + b) / 2;val = f(alpha);end步骤4:实现Powell法迭代算法现在,我们可以基于上述步骤定义Powell法的主迭代算法。
机械优化设计程序
机械优化设计程序(2009-04-24 19:30:52)转载标签:机械优化设计程序powell外点惩罚函数法c杂谈今天终于交了机械优化设计作业,程序贴出来,那位要用就拿去吧,资源共享了,O(∩_∩)O下面是利用外点惩罚函数(Powell)法求解三维目标函数最优解与最优值的程序,用C++编写。
#include<iostream.h>#include<math.h>double lamta[10]={0, 1.0 ,0 ,0 ,0 ,1 ,0 ,0 ,0 ,1};//鲍威尔方法初始化方向,线性无关double lamta1[3]={0, 0 , 0};//暂存新的搜索方向double x1[4]={0, 0 ,0, 0 };//x1到x3用于存储各共轭方向的点double x2[4]={0, 0 ,0, 0 };double x3[4]={0, 0 ,0, 0 };double x4[4]={0, 0 ,0, 0 };//x4用于中间判断double x5[4]={0, 0 ,0, 0 };//x5用存放于更换方向后产生的新点int m=0;//标志double x_[4]={0, 0, 0, 0};//暂存鲍威尔最优解double x0[4]={0, 2, 2 , 2};//初值double c=10;//递减系数double e=0.00001;//精度控制double r0=1;//初始惩罚因子double r=1;//函数声明部分void Powell(double r); //鲍威尔方法函数double fxy(double x1,double x2,double x3,double r); //待求函数double ysearch(double x); //一维搜索的目标函数void search(double &a,double &b,double h); //区间搜索double yellowcut(double &a,double &b); //黄金分割void sort(double *p,int size);//选择法排序void main() //约束优化方法主函数入口{cout<<"请输入精度"<<endl;cin>>e;changyan:Powell(r);double cmpare[4];int flag1=0;for (int i=1;i<=3;i++){cmpare[i]=x_[i]-x0[i];if (fabs(cmpare[i])<e){flag1++;}}if (flag1==3){cout<<"最优解为:"<<"x1="<<x_[1]<<" "<<"x2="<<x_[2]<<" "<<"x3="<<x_[3]<<endl;cout<<"最小值为"<<fxy(x_[1],x_[2],x_[3],r)<<endl;}else{for (int j=1;j<=3;j++){x0[j]=x_[j];}r=c*r;goto changyan;}}//子函数定义部分double fxy(double x1,double x2,double x3,double r)//待求函数{double m,n,p;m=(-x1>0)?(-x1):0;n=(-x2>0)?(-x2):0;p=(-x3>0)?(-x3):0;return //惩罚函数1000-x1*x1-2*x2*x2-x3*x3-x1*x2-x1*x3+r*(m*m+n*n+p*p)+r*((x1*x1+x2*x2+x3*x3-25)*(x1*x1+x2*x2+x3 *x3-25)+(8*x1+14*x2+7*x3-56)*(8*x1+14*x2+7*x3-56));}void Powell(double r) //鲍威尔方法函数定义double det=0.0001; //迭代精度int k;my1: for (k=1;k<=3;k++){m=3*k-2;double a=0,b=0,xo=0;search(a,b,1); //完成区间搜索double temp;temp=yellowcut(a,b);//黄金分割法int n=3*k-2;for (int i=1;i<=3;i++){switch (k){case 1:x1[i]=x0[i]+temp*lamta[n++];break;case 2:x2[i]=x1[i]+temp*lamta[n++];break;case 3:x3[i]=x2[i]+temp*lamta[n++];break;default :break;}}}double cmp[4];int flag=0;for (int i=1;i<=3;i++){cmp[i]=x3[i]-x0[i];if (fabs(cmp[i])<det){flag++;}}if (flag==3) //找到最优解x_[1]=x3[1];x_[2]=x3[2];x_[3]=x3[3];}else{double fy[4];fy[0]=fxy(x0[1],x0[2],x0[3],r);fy[1]=fxy(x1[1],x1[2],x1[3],r);fy[2]=fxy(x2[1],x2[2],x2[3],r);fy[3]=fxy(x3[1],x3[2],x3[3],r); double fyy[3]; for (int ii=0;ii<3;ii++){fyy[ii]=fy[ii]-fy[ii+1];}sort(fyy,3);for (ii=1;ii<=3;ii++){x4[ii]=2*x3[ii]-x0[ii];}double f0,f3,f4;f0=fy[0];f3=fy[3];f4=fxy(x4[1],x4[2],x4[3],r);if ((f0+f4-2*f3)/2>=fyy[2]){if (f3<f4){for (int t=1;t<=3;t++){x0[t]=x3[t];}}elsefor (int t=1;t<=3;t++){x0[t]=x4[t];}}goto my1;}else{for (int t=0;t<3;t++){lamta1[t]=x3[t+1]-x0[t+1];}m=0; //switch 标志!double aa=0,bb=0;search(aa,bb,1);double temp1;temp1=yellowcut(aa,bb);for (int i=1;i<=3;i++){x5[i]=x3[i]+temp1*lamta1[i-1];}for (i=1;i<=3;i++){x0[i]=x5[i];}for (i=1;i<=6;i++){lamta[i]=lamta[i+3];}for (i=1;i<=3;i++){lamta[6+i]=lamta1[i-1];}goto my1;}}}double ysearch(double x) //一维搜索的目标函数{switch (m){case 1: return fxy(x0[1]+x*lamta[m],x0[2]+x*lamta[m+1],x0[3]+x*lamta[m+2],r);break;case 4: return fxy(x1[1]+x*lamta[m],x1[2]+x*lamta[m+1],x1[3]+x*lamta[m+2],r);break;case 7: return fxy(x2[1]+x*lamta[m],x2[2]+x*lamta[m+1],x2[3]+x*lamta[m+2],r);break;case 0: return fxy(x3[1]+x*lamta1[0],x3[2]+x*lamta1[1],x3[3]+x*lamta1[2],r);break;//更改方向后的一维搜索default:return 0; break;}}void search(double &a,double &b,double h) //区间搜索{double a1,a2,a3,y1,y2,y3;h=1;a1=a,y1=ysearch(a1);a2=a+h,y2=ysearch(a2);if(y2>=y1){h=-h,a3=a1,y3=y1;a1=a2,y1=y2,a2=a3,y2=y3;}a3=a2+h,y3=ysearch(a3);while(y3<=y2){h=2*h;a1=a2,y1=y2,a2=a3,y2=y3;a3=a2+h,y3=ysearch(a3);}if(h<0)a=a3,b=a1;else a=a1,b=a3;}double yellowcut(double &a,double &b){double e; //黄金分割法求解e=0.001;double c,fc;c=a+0.382*(b-a);fc=ysearch(c);double d,fd;double xo;d=a+0.618*(b-a);fd=ysearch(d);label2: if (fc<=fd){b=d;d=c;fd=fc;c=a+0.382*(b-a);fc=ysearch(c);}else{a=c;c=d;fc=fd;d=a+0.618*(b-a);fd=ysearch(d);}if ((b-a)<=e){xo=(a+b)/2;}elsegoto label2;return xo;}void sort(double *p,int size){//选择法排序int i,j;double k;for(i=0;i<size-1;i++)for(j=i+1;j<size;j++)if(*(p+i)>*(p+j)){k=*(p+i);*(p+i)=*(p+j);*(p+j)=k;} }。
机械优化设计鲍威尔法
机械优化设计鲍威尔法
机械优化设计鲍威尔法(Powell method)是一种常用的非线性优化
算法,它以鲍威尔的名字命名,用于解决无约束非线性优化问题。
该方法
在各个领域都有广泛的应用,如工程优化设计、机器学习等。
下面将详细
介绍机械优化设计鲍威尔法的原理及应用。
鲍威尔法的具体步骤如下:
1.初始化参数:选择初始设计参数和方向。
2.寻找一维极小值点:沿着方向找到目标函数在该方向上的极小值点。
3.更新方向:通过比较前后两个极小值点的差异来更新方向。
4.迭代优化:重复步骤2和步骤3,直到达到指定的收敛条件。
鲍威尔法的优点是收敛速度较快、计算量较小,同时可以处理非线性
的优化问题。
然而,该方法也存在一些不足之处,如可能陷入局部最优解、对初值敏感等。
机械优化设计鲍威尔法在工程领域中有广泛的应用。
例如,在机械结
构设计中,可以利用鲍威尔法来优化结构参数,以满足特定的性能指标。
在汽车工业中,可以使用鲍威尔法来优化车辆的燃油效率和性能。
在航空
航天领域,可以利用该方法来优化飞行器的飞行性能。
此外,该方法还可
以用于机器学习中的参数优化,如调整神经网络的权重和偏置等。
总之,机械优化设计鲍威尔法是一种常用的非线性优化算法,通过迭
代逼近最优解。
虽然该方法有一些不足之处,但在实际应用中具有广泛的
适用性,尤其在工程优化设计和机器学习等领域。
通过使用该方法,可以
优化设计参数,改进性能指标,提高工程效率和产品质量。
机械优化设计鲍威尔法编程
}
}
// 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)法是直接利用函数值 来构造共扼方向的一种共扼方向法。基本 思想是在不用求导数的前提下,在迭代中 逐次构造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个线性无关的向量
0.618法牛顿法Powell法优化设计MATLAB程序设计
机械优化设计*名:***学号:************ 班级:机械工程15班1.用0.618法求一元函数f(t)=t^2+2*t在区间[-3,5]中的极小点,要求计算到区间的长度小于0.05。
用MATLAB编程如下:function [x,minf]=minHJ(f,a,b,eps)format long;if nargin==3eps=0.05;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);if fl>fua=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;x=NaN;minf=NaN;return;endx=(a+b)/2;minf=subs(f,findsym(f),x);minf = double(minf);format short;syms t;f=t^2+2*t;[x,fx]=minHJ(f,-3,5,eps)运行结果:2.用牛顿法求函数f = x1^2 + 4*x2^2 + x3^2 - 2*x1 + 18*x3的极小点。
用MATLAB编程如下:function [diff_x] = function_c1(x1,x2,x3)syms x1x2x3f = x1^2 + 4*x2^2 + x3^2 - 2*x1 + 18*x3;diff_x1 = diff(f,x1);diff_x2 = diff(f,x2);diff_x3 = diff(f,x3);diff_x = [diff_x1; diff_x2; diff_x3];returnendfunction [H] = function_c2(x1,x2,x3)syms x1x2x3f = x1^2 + 4*x2^2 + x3^2 - 2*x1 + 18*x3;diff_x1 = diff(f,x1);diff_x2 = diff(f,x2);diff_x3 = diff(f,x3);diff_x1_x1 = diff(diff_x1,x1);diff_x2_x2 = diff(diff_x2,x2);diff_x3_x3 = diff(diff_x3,x3);diff_x1_x2 = diff(diff_x1,x2);diff_x1_x3 = diff(diff_x1,x3);diff_x2_x1 = diff_x1_x2;diff_x2_x3 = diff(diff_x2,x3);diff_x3_x1 = diff_x1_x3;diff_x3_x2 = diff_x2_x3;H = [diff_x1_x1 diff_x1_x2 diff_x1_x3;...diff_x2_x1 diff_x2_x2 diff_x2_x3;...diff_x3_x1 diff_x3_x2 diff_x3_x3];returnende = 1;x1 = 0;x2 = 0;x3 = 0;while e>0.00001f_d = function_c1(x1,x2,x3);f = subs(f_d, 'x1,x2,x3', [x1,x2,x3]);H = function_c2(x1,x2,x3);H_inv = inv(H);x = [x1;x2;x3] - H_inv * f;e = abs(x(1)-x1) / x(1);x1 = x(1);x2 = x(2);x3 = x(3);enddisp(x)运行结果:3.用Powell法求函数f=t^2+s^2-t*s-10*t-4*s+60的最优点X*=[t,s]T。
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.在当前方向上进行线,找到使目标函数值下降的步长。
4.更新设计参数,将方向和步长相乘加到当前设计参数上。
5.检查当前设计参数的变化是否满足终止条件,如果满足则结束,否则返回步骤26.返回最佳设计参数组合及对应的目标函数值。
下面以一个简单的机械优化设计案例为例进行详细说明。
假设有一个弹簧悬挂系统,需要设计合适的弹簧刚度和阻尼系数来满足特定的振动要求。
目标函数为最小化系统振动幅值。
首先,需要定义设计参数和目标函数。
设计参数可以选择弹簧刚度和阻尼系数。
目标函数可以定义为系统振动幅值的平方。
其次,需要确定方向和初始步长。
对于弹簧刚度和阻尼系数来说,方向可以选择正方向和负方向。
初始步长可以根据经验或试验来确定。
然后,根据鲍威尔法编程的步骤,进行迭代。
首先,初始化设计参数和步长。
然后,计算当前设计参数下的目标函数值。
接下来,根据当前方向和步长进行线,找到使目标函数值下降的步长。
再更新设计参数,将方向和步长相乘加到当前设计参数上。
最后,检查当前设计参数的变化是否满足终止条件,如果满足则结束,否则返回到线步骤继续。
最后,得到最佳的设计参数组合及对应的目标函数值。
可以根据实际情况进行设计参数的调整和优化,以获得更好的机械性能。
总之,鲍威尔法编程是一种常用的机械优化设计方法。
通过不断迭代的方式,寻找最佳的设计参数组合,以满足设计要求。
其原理和步骤相对简单,但需要结合具体问题进行参数的选择和调整。
机械优化设计-鲍威尔法
机械优化设计-鲍威尔法#include <stdio.h>#include <math.h>#define m 10 /*数组长度m >= 维数n */ float f(float x[]);void mjtf(int n,float x0[],float h,float s[],float a[],float b[]);void mhjfgf(int n,float a[],float b[],float flag,float x[]);void mbwef(int n,float x0[],float h,float flag,float a[],floatb[],float x[]);/*目标函数(n维)*//*入口参数:x :n维数组,自变量 *//*返回值 :函数值 */float f(float x[]){float result;result=60-10*x[0]-4*x[1]+x[0]*x[0]+x[1]*x[1]-x[0]*x[1];//result=10*(float)pow((x[0]+x[1]-5),2)+(float)pow((x[0]-x[1]),2);return result;}/*多维进退法子程序*//*入口参数:n :优化模型维数x0 :n维数组,初始点坐标h :初始搜索步长s :n维数组,搜索方向 *//*出口参数:a :n维数组,搜索区间下限b :n维数组,搜索区间上限*/void mjtf(int n,float x0[],float h,float s[],float a[],float b[]) {int i;float x1[m],x2[m],x3[m],f1,f2,f3;for(i=0;i<n;i++) /*计算初始两试点*/{x1[i]=x0[i];x2[i]=x0[i]+h*s[i];}f1=f(x1);f2=f(x2);if(f2>=f1) /*判断搜索方向*/{ /*搜索方向为反向,转身*/h=(-1)*h;for(i=0;i<n;i++)x3[i]=x1[i];f3=f1;for(i=0;i<n;i++)x1[i]=x2[i];f1=f2;for(i=0;i<n;i++)x2[i]=x3[i];f2=f3;} /*搜索方向为正向*/for(i=0;i<n;i++) /*计算第三试点*/x3[i]=x2[i]+h*s[i];f3=f(x3);while(f3<f2) /*判断是否未完成搜索*/ { /*未完成,继续搜索*/h=2*h;for(i=0;i<n;i++)x1[i]=x2[i];f1=f2;for(i=0;i<n;i++)x2[i]=x3[i];f2=f3;for(i=0;i<n;i++)x3[i]=x2[i]+h*s[i];f3=f(x3);} /*已完成*/for(i=0;i<n;i++) /*输出初始搜索区间*/ {if(h>0)// if(x1[i]<x3[i]){a[i]=x1[i];b[i]=x3[i];}else{a[i]=x3[i];b[i]=x1[i];}}}/*多维黄金分割法子程序*//*入口参数:n :优化模型维数a :n维数组,搜索区间下限b :n维数组,搜索区间上限flag :迭代精度 */ /*出口参数:x :n维数组,极小点坐标 */void mhjfgf(int n,float a[],float b[],float flag,float x[]) {int i;float x1[m],x2[m],f1,f2,sum;for(i=0;i<n;i++) /*计算初始两试点*/x1[i]=b[i]-(float)0.618*(b[i]-a[i]); f1=f(x1);for(i=0;i<n;i++)x2[i]=a[i]+(float)0.618*(b[i]-a[i]); f2=f(x2);do{if(f1<=f2) /*判断消去区间*/{ /*消去右*/for(i=0;i<n;i++)b[i]=x2[i];for(i=0;i<n;i++)x2[i]=x1[i];f2=f1;for(i=0;i<n;i++)x1[i]=b[i]-(float)0.618*(b[i]-a[i]); f1=f(x1);}else{ /*消去左*/for(i=0;i<n;i++)a[i]=x1[i];for(i=0;i<n;i++)x1[i]=x2[i];f1=f2;for(i=0;i<n;i++)x2[i]=a[i]+(float)0.618*(b[i]-a[i]);f2=f(x2);}sum=0;for(i=0;i<n;i++)sum+=(b[i]-a[i])*(b[i]-a[i]);}while(sqrt(sum)>flag*0.1); /*判断是否未达到精度要求,若未达到则返回继续缩小区间*/for(i=0;i<n;i++)x[i]=(float)0.5*(b[i]+a[i]); /*已达到,输出极小点*/ }/*鲍威尔法子程序*//*入口参数:n :优化模型维数x0 :n维数组,初始点坐标h :初始搜索步长flag :迭代精度a :n维数组,搜索区间下限b :n维数组,搜索区间上限*//*出口参数:x :n维数组,极小点坐标 */void mbwef(int n,float x0[],float h,float flag,float a[],floatb[],float x[]){int i,j,k,r;float x1[m],x2[m],f0,f1,f2,fn[m],s[m][m],sum;for(i=0;i<n;i++) /*方向矩阵初始化*/for(k=0;k<n;k++)if(i==k)s[i][k]=1;elses[i][k]=0;k=1;while(1){for(i=0;i<n;i++)x1[i]=x0[i];for(i=0;i<n;i++) /*依次按每个方向搜索*/{mjtf(n,x1,h,s[i],a,b); /*调用多维进退法子程序*/mhjfgf(n,a,b,flag,x1); /*调用多维黄金分割法子程序*/fn[i]=f(x0)-f(x1); /*计算函数下降值*/}for(i=0;i<n;i++) /*计算映射点*/x2[i]=2*x1[i]-x0[i];for(i=1;i<n;i++) /*找出函数下降值中的最大值存入fn[0]*/ if(fn[0]<fn[i]){fn[0]=fn[i];r=i;}elser=0;f0=f(x0); /*计算始点、终点和映射点的函数值*/f1=f(x1);f2=f(x2);if(f2>=f0||(f0-2*f1+f2)*(f0-f1-fn[0])*(f0-f1-fn[0])>=0.5*fn[0]*(f0-f2)*(f0-f2)){ /*判断是否需要换方向*//*不需要换*/sum=0; /*计算一轮中终点与始点的距离*/for(i=0;i<n;i++)sum+=(x1[i]-x0[i])*(x1[i]-x0[i]);if(f1<=f2) /*判断将终点还是将映射点赋给下一轮始点*/for(i=1;i<n;i++) /*将终点赋给下一轮始点*/x0[i]=x1[i];elsefor(i=1;i<n;i++) /*将映射点赋给下一轮始点*/x0[i]=x2[i];}else{ /*需要换*/for(i=r;i<n;i++) /*去掉第r+1个方向,其它方向前移*/ for(j=0;j<n;j++)s[i][j]=s[i+1][j];for(i=0;i<n;i++) /*生成新方向放在最后*/s[n][i]=x1[i]-x0[i];mjtf(n,x1,h,s[n],a,b); /*按新方向搜索一次*/mhjfgf(n,a,b,flag,x1);sum=0; /*计算一轮中终点与始点的距离*/for(i=0;i<n;i++)sum+=(x1[i]-x0[i])*(x1[i]-x0[i]);for(i=0;i<n;i++) /*将终点赋给下一轮始点*/x0[i]=x1[i];}if(sqrt(sum)<=flag) /*判断是否满足精度要求*/ break;elsek+=1;}for(i=0;i<n;i++) /*输出极小点坐标*/x[i]=x1[i];}/*鲍威尔法主程序*/void main(){int i,n;float h,flag,x0[m],a[m],b[m],x[m]; printf("\n<鲍威尔法>\n");printf("请输入维数:\n");scanf("%d",&n);printf("请输入初始点:");for(i=0;i<n;i++){printf("\nx0[%d]=",i);scanf("%f",&x0[i]);}printf("\n请输入初始步长:\n"); scanf("%f",&h);printf("\n请输入精度:\n");scanf("%f",&flag);mbwef(n,x0,h,flag,a,b,x);printf("\n极小点坐标为:\n");for(i=0;i<n;i++)printf("x[%d]=%f\n",i,x[i]);printf("\n极小值为:\n%f\n",f(x)); }。
用powell法优化设计程序与一维搜索黄金分割法组合
用Powell法优化设计程序与一维搜索黄金分割法组合编程求解函数f(x) x1 2x2 4x-| 2x-|X2的极小点X,初始点X0=[1,1]T,迭代精度F0.001o 解:已知f(x) x2 2x| 4x-| 2x1x2,初始点x°=[1,1]T,迭代精度卩0.001。
在该优化设计过程中,黄金分割搜索法作为POWELL算法主程序中的一部分。
在POWELL算法运行过程中会多次调用黄金分割搜索算法程序。
这样可以缩短优化设计计算时间。
1. MATLAB源程序代码1.1关于a的目标函数源代码1.2 一元函数最小值区间函数源代码endif alpha1>alpha3tem = alpha1;alpha1 = alpha3;alpha3 = tem;a = alpha1;b = alpha3;else a=alpha1;b = alpha3;end1.3黄金分割搜索法函数源代码1.4 POWELL算法程序源代码2. 运行程序计算>> f = @(x1,x2) x1A2+2*x2A2-4*x1-2*x1*x2f =@(x1,x2)x1A2+2*x2A2-4*x1-2*x1*x2>> [z,fmi n] = powell(f)z =3.99981.9998fmin =-8.0000>>经过计算可知,极小值点为(3.9998,1.9998),极小值为-83. 验算通过MATLAB软件内置fminsearch函数进行验算。
3.1绘图通过观察可以发现极小值点在(4,2)附近3.2 fmi nsearch 函数验算>> f = @(x) x⑴A2+2*x (2) A2-4*x(1)-2*x(1)*x(2)f =@(x)x(1)A2+2*x (2) A2-4*x(1)-2*x(1)*x(2)>> [x,fval] = fmi nsearch(f,[1,1])x =4.0000 2.0000fval =-8.0000>>通过fminsearch函数验算,确认极小值点为(4,2),极小值为-8。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
#include "stdio.h"
#include "stdlib.h"
#include "math.h"
double objf(double x[])
{double ff;
ff=x[0]*x[0]+x[1]*x[1]-x[0]*x[1]-10*x[0]-4*x[1]+60;
return(ff);
}
double gold(double a[],double b[],double eps,int n,double xx[]) {int i;
double f1,f2,*x[2],ff,q,w;
for(i=0;i<2;i++)
x[i]=(double *)malloc(n*sizeof(double));
for(i=0;i<n;i++)
{*(x[0]+i)=a[i]+0.618*(b[i]-a[i]);
*(x[1]+i)=a[i]+0.382*(b[i]-a[i]);
}
f1=objf(x[0]);
f2=objf(x[1]);
do
{if(f1>f2)
{for(i=0;i<n;i++)
{b[i]=*(x[0]+i);
*(x[0]+i)=*(x[1]+i);
}
f1=f2;
for(i=0;i<n;i++)
*(x[1]+i)=a[i]+0.382*(b[i]-a[i]);
f2=objf(x[1]);
}
else
{ for(i=0;i<n;i++)
{a[i]=*(x[1]+i);
*(x[1]+i)=*(x[0]+i);}
f2=f1;
for(i=0;i<n;i++)
*(x[0]+i)=a[i]+0.618*(b[i]-a[i]);
f1=objf(x[0]);
}
q=0;
for(i=0;i<n;i++)
q=q+(b[i]-a[i])*(b[i]-a[i]);
}while(w>eps);
for(i=0;i<n;i++)
xx[i]=0.5*(a[i]+b[i]);
ff=objf(xx);
for(i=0;i<2;i++)
free(x[i]);
return(ff);
}
void jtf(double x0[],double h0,double s[],int n,double a[],double b[]) {int i;
double *x[3],h,f1,f2,f3;
for(i=0;i<3;i++)
x[i]=(double *)malloc(n*sizeof(double));
h=h0;
for(i=0;i<n;i++)
*(x[0]+i)=x0[i];
f1=objf(x[0]);
for(i=0;i<n;i++)
*(x[1]+i)=*(x[0]+i)+h*s[i];
f2=objf(x[1]);
if(f2>=f1)
{h=-h0;
for(i=0;i<n;i++)
*(x[2]+i)=*(x[0]+i);
f3=f1;
for(i=0;i<n;i++)
{*(x[0]+i)=*(x[1]+i);
*(x[1]+i)=*(x[2]+i);
}
f1=f2;
f2=f3;
}
for(;;)
{h=2*h;
for(i=0;i<n;i++)
*(x[2]+i)=*(x[1]+i)+h*s[i];
f3=objf(x[2]);
if(f2<f3) break;
else
{ for(i=0;i<n;i++)
{*(x[0]+i)=*(x[1]+i);
*(x[1]+i)=*(x[2]+i);
f1=f2;
f2=f3;
}
}
if(h<0)
for(i=0;i<n;i++)
{a[i]=*(x[2]+i);
b[i]=*(x[0]+i);
}
else
for(i=0;i<n;i++)
{a[i]=*(x[0]+i);
b[i]=*(x[2]+i);
}
for(i=0;i<3;i++)
free(x[i]);
}
double oneoptim(double x0[],double s[],double h0,double epsg,int n,double x[]) {double *a,*b,ff;
a=(double *)malloc(n*sizeof(double));
b=(double *)malloc(n*sizeof(double));
jtf(x0,h0,s,n,a,b);
ff=gold(a,b,epsg,n,x);
free(a);
free(b);
return (ff);
}
double powell(double p[],double h0,double eps,double epsg,int n,double x[]) {int i,j,m;
double *xx[4],*ss,*s;
double f,f0,f1,f2,f3,fx,dlt,df,sdx,q,d;
ss=(double *)malloc(n*(n+1)*sizeof(double));
s=(double *)malloc(n*sizeof(double));
for(i=0;i<n;i++)
{for(j=0;j<=n;j++)
*(ss+i*(n+1)+j)=0;
*(ss+i*(n+1)+i)=1;
}
for(i=0;i<4;i++)
xx[i]=(double *)malloc(n*sizeof(double));
for(i=0;i<n;i++)
*(xx[0]+i)=p[i];
for(;;)
{for(i=0;i<n;i++)
{*(xx[1]+i)=*(xx[0]+i);
x[i]=*(xx[1]+i);
}
f0=f1=objf(x);
dlt=-1;
for(j=0;j<n;j++)
{for(i=0;i<n;i++)
{*(xx[0]+i)=x[i];
*(s+i)=*(ss+i*(n+1)+j);
}
f=oneoptim(xx[0],s,h0,epsg,n,x);
df=f0-f;
if(df>dlt)
{dlt=df;
m=j;
}
}
sdx=0;
for(i=0;i<n;i++)
sdx=sdx+fabs(x[i]-(*(xx[1]+i)));
if(sdx<eps)
{free(ss);
free(s);
for(i=0;i<4;i++)
free(xx[i]);
return(f);
}
for(i=0;i<n;i++)
*(xx[2]+i)=x[i];
f2=f;
for(i=0;i<n;i++)
{*(xx[3]+i)=2*(*(xx[2]+i)-(*(xx[1]+i)));
x[i]=*(xx[3]+i);
}
fx=objf(x);
f3=fx;
q=(f1-2*f2+f3)*(f1-f2-dlt)*(f1-f2-dlt);
d=0.5*dlt*(f1-f3)*(f1-f3);
if((f3<f1)||(q<d))
{if(f2<=f3)
for(i=0;i<n;i++)
*(xx[0]+i)=*(xx[2]+i);
else
for(i=0;i<n;i++)
*(xx[0]+i)=*(xx[3]+i);
}
else
{for(i=0;i<n;i++)
{*(ss+(i+1)*(n+1))=x[i]-(*(xx[1]+i));
*(s+i)=*(ss+(i+1)*(n+1));
}
f=oneoptim(xx[0],s,h0,epsg,n,x);
for(i=0;i<n;i++)
*(xx[0]+i)=x[i];
for(j=m+1;j<=n;j++)
for(i=0;i<n;i++)
*(ss+i*(n+1)+j-1)=*(ss+i*(n+1)+j);
}
}
}
void main()
{double p[]={1,2};
double ff,x[2];
ff=powell(p,0.3,0.001,0.0001,2,x);
printf("x[0]=%f,x[1]=%f,ff=%f\n",x[0],x[1],ff); getchar();
}。