鲍威尔法C++源程序
- 格式:doc
- 大小:42.50 KB
- 文档页数:6
#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]);w=sqrt(q);}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(;;){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);elsefor(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);}}}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); getch();}。
实验报告实验名称:鲍威尔法院(系):机电学院专业班级:机械制造及其自动化姓名:学号: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的基本操作方法,了解鲍威尔法的原理与基本运用。
机械优化设计鲍威尔法编程鲍威尔法(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)为示例函数,可以根据具体的优化设计问题进行修改。
include <>include <>include <>include <>include <>include <>define n1 2define ttdefine ad//定义常量//tt 初始迭代步长//ad 收敛精度float ia;float fnyfloat x{float f;f=10powx0+x1-5,2+powx0-x1,2; //目标函数returnf;}float iteratefloat x,float a, float s{float x1;x1=float mallocn1 sizeoffloat;for int i=0;i<n1;i++x1i=xi+asi;returnx1;}float funcfloat x,float a,float s{float x1;x1=iteratex,a,s;float f=fnyx1;returnf;}void findingfloat a3,float f3,float xk,float s {float t=tt;float a1,f1;a0=0; f0=funcxk,a0,s;for int i=0;;i++{a1=a0+t; f1=funcxk,a1,s;if f1<f0break;if fabsf1-f0>=ad{t=-t;a0=a1; f0=f1;}else{if ia==1 return;t=t/2; ia=1;}}for i=0;;i++{a2=a1+t; f2=funcxk,a2,s;if f2>f1break;t=2t;a0=a1; f0=f1;a1=a2; f1=f2;}if a0>a2{a1=a0; f1=f0;a0=a2; f0=f2;a2=a1; f2=f1;}return;}//second insertfloat lagrangefloat xk,float ft,float s{float a3,f3;float b,c,d,aa;finding a,f,xk,s;for int i=0;;i++{if ia==1{aa=a1; ft=f1;break;}d=powa0,2-powa2,2a0-a1-powa0,2-powa1,2a0-a2;iffabsd==0break;c=f0-f2a0-a1-f0-f1a0-a2/d;iffabsc==0break;b=f0-f1-cpowa0,2-powa1,2/a0-a1;aa=-b/2c;ft=funcxk,aa,s;if fabsaa-a1<=ad{if ft>f1 aa=a1;break;}if aa>a1{if ft>f1{a2=aa; f2=ft;}else if ft<f1{a0=a1; a1=aa;f0=f1; f1=ft;}else if ft==f1{a2=aa; a0=a1;f2=ft; f0=f1;a1=a0+a2/2; f1=funcxk,a1,s;}}else{if ft>f1{a0=aa; f0=ft;}else if ft<f1{a2=a1; a1=aa;f2=f1; f1=ft;}else if ft==f1{a0=aa; a2=a1;f0=ft; f2=f1;a1=a0+a2/2; f1=funcxk,a1,s;}}}if ft>f1{aa=a1; ft=f1;}return aa;}float powellfloat xk{float hn1n1,sn1={0,0},ffn1+1={0,0,0};float f1,f3,aa;float dkn1,x0,xk1n1;int m=0,i,j;for i=0;i<n1;i++{forj=0;j<n1;j++{hij=0;if j==ihij=1;}}for int k=0;;k++{ff0=fnyxk;x0=xk;for i=0;i<n1;i++{for j=0;j<n1;j++sj=hij;float aa=lagrangexk,&ffi+1,s;xk=iteratexk,aa,s;}for i=0;i<n1;i++{float a,b;dki=xki-x0i;xk1i=2xki-x0i;}float max=fabsff1-ff0;for i=1;i<n1;i++if fabsffi+1-ffi>max{max=fabsffi+1-ffi;m=i;}f3=fnyxk1;if f3<ff0&&ff0+f3-2ff2powff0-ffn1-max,2<maxpowff0-f3,2 {aa=lagrangexk,&f1,dk;xk=iteratexk,aa,dk;for i=m;i<n1-1;i++for j=0;j<n1;j++hij=hi+1j;for j=0;j<n1;j++hn1-1j=dkj;}else{ifffn1>=f3 xk=xk1;}float xq=0;for i=0;i<n1;i++xq+=powxki-x0i,2;if xq<=adbreak;}returnxk;}void main{float xkn1={0,0};//取初始点float xx;xx=float mallocn1 sizeoffloat;xx=powellxk;float ff=fnyxx;cout<<"优化的结果为:"<<endl;printf"\n\nThe Optimal Design Result Is:\n";for int i=0;i<n1;i++printf"\n\t x%d =%f",i+1,xxi;printf"\n\t f=%f",ff;getch;}。
10种简单的数字滤波C语言源程序算法(2009-11-09 10:25:08)假定从8位AD中读取数据(如果是更高位的AD可定义数据类型为int),子程序为get_ad();1、限幅滤波法(又称程序判断滤波法)A、方法:根据经验判断,确定两次采样允许的最大偏差值(设为A)每次检测到新值时判断:如果本次值与上次值之差<=A,则本次值有效如果本次值与上次值之差>A,则本次值无效,放弃本次值,用上次值代替本次值B、优点:能有效克服因偶然因素引起的脉冲干扰C、缺点无法抑制那种周期性的干扰平滑度差2、中位值滤波法A、方法:连续采样N次(N取奇数)把N次采样值按大小排列取中间值为本次有效值B、优点:能有效克服因偶然因素引起的波动干扰对温度、液位的变化缓慢的被测参数有良好的滤波效果C、缺点:对流量、速度等快速变化的参数不宜3、算术平均滤波法A、方法:连续取N个采样值进行算术平均运算N值较大时:信号平滑度较高,但灵敏度较低N值较小时:信号平滑度较低,但灵敏度较高N值的选取:一般流量,N=12;压力:N=4B、优点:适用于对一般具有随机干扰的信号进行滤波这样信号的特点是有一个平均值,信号在某一数值范围附近上下波动C、缺点:对于测量速度较慢或要求数据计算速度较快的实时控制不适用比较浪费RAM4、递推平均滤波法(又称滑动平均滤波法)A、方法:把连续取N个采样值看成一个队列队列的长度固定为N每次采样到一个新数据放入队尾,并扔掉原来队首的一次数据.(先进先出原则) 把队列中的N个数据进行算术平均运算,就可获得新的滤波结果N值的选取:流量,N=12;压力:N=4;液面,N=4~12;温度,N=1~4 B、优点:对周期性干扰有良好的抑制作用,平滑度高适用于高频振荡的系统C、缺点:灵敏度低对偶然出现的脉冲性干扰的抑制作用较差不易消除由于脉冲干扰所引起的采样值偏差不适用于脉冲干扰比较严重的场合比较浪费RAM5、中位值平均滤波法(又称防脉冲干扰平均滤波法)A、方法:相当于“中位值滤波法”+“算术平均滤波法”连续采样N个数据,去掉一个最大值和一个最小值然后计算N-2个数据的算术平均值N值的选取:3~14B、优点:融合了两种滤波法的优点对于偶然出现的脉冲性干扰,可消除由于脉冲干扰所引起的采样值偏差C、缺点:测量速度较慢,和算术平均滤波法一样比较浪费RAM6、限幅平均滤波法A、方法:相当于“限幅滤波法”+“递推平均滤波法”每次采样到的新数据先进行限幅处理,再送入队列进行递推平均滤波处理B、优点:融合了两种滤波法的优点对于偶然出现的脉冲性干扰,可消除由于脉冲干扰所引起的采样值偏差C、缺点:比较浪费RAM7、一阶滞后滤波法A、方法:取a=0~1本次滤波结果=(1-a)*本次采样值+a*上次滤波结果B、优点:对周期性干扰具有良好的抑制作用适用于波动频率较高的场合C、缺点:相位滞后,灵敏度低滞后程度取决于a值大小不能消除滤波频率高于采样频率的1/2的干扰信号8、加权递推平均滤波法A、方法:是对递推平均滤波法的改进,即不同时刻的数据加以不同的权通常是,越接近现时刻的数据,权取得越大。
机械优化设计——鲍威尔法机械优化设计班级:0841001成员:张波2010213217张建2010213214潘阳瑞20102132272013年6月鲍威尔法鲍威尔(Powell)法是直接利用函数值来构造共轭方向的一种方法。
基本思想:在不用导数的前提下,在迭代中逐次构造G 的共轭方向。
一(基本算法:(二维情况描述鲍威尔的基本算法)0T1)任选一初始点x,再选两个线性无关的向量,如坐标轴单位向量e=[1,0]和1T=[0,1]作为初始搜索方向。
e20002)从出发,顺次沿、作一维搜索,得、点,两点连线得一新 xeexx12121001方向 d,x,xd2011 用代替e形成两个线性无关向量,e,作为下一轮迭代的搜索方向。
再从xdd1,1201出发,沿作一维搜索得点,作为下一轮迭代的初始点。
xd111113)从出发,顺次沿、作一维搜索,得到点、,两点连线得一新方向: exxxd122211。
d,x,x21*22沿作一维搜索得点,即是二维问题的极小点。
xdx把二维情况的基本算法扩展到n维,则鲍威尔基本算法的要点是:在每一轮迭代中总有一个始点(第一轮的始点是任选的初始点)和n个线性独立的搜索方向。
从始点出发顺次沿n个方向作一维搜索得一终点,由始点和终点决定了一个新的搜索方向。
用这个方向替换原来n个方向中的一个,于是形成新的搜索方向组。
替换的原则是去掉原方向组的第一个方向而将新方向排在原方向的最后。
此外规定,从这一轮的搜索终点出发沿新的搜索方向作一维搜索而得到的极小点,作为下一轮迭代的始点。
这样就形成算法的循环。
图1.二维情况下的鲍威尔法二(改进算法在鲍威尔基本算法中,每一轮迭代都用连结始点和终点所产生出的搜索方向去替换原向量组中的第一个向量,而不管它的“好坏”,这是产生向量组线性相关的原因所在。
在改进的算法中首先判断原向量组是否需要替换。
如果需要替换,还要进一步判断原向量组中哪个向量最坏,然后再用新产生的向量替换这个最坏的向量,以保证逐次生成共轭方向。
C语言经典源程序100例1. Hello, World!这是C语言中最基本的程序,用于显示"Hello, World!"。
```c#include <stdio.h>int main() {printf("Hello, World!\n");return 0;}```2. 计算两数之和这个程序用于计算两个整数的和,并将结果输出。
```c#include <stdio.h>int main() {int num1, num2, sum;printf("请输入两个整数:");scanf("%d %d", &num1, &num2);sum = num1 + num2;printf("两数之和为:%d\n", sum);return 0;}```3. 判断奇偶数这个程序用于判断一个整数是奇数还是偶数。
```c#include <stdio.h>int main() {int num;printf("请输入一个整数:");scanf("%d", &num);if (num % 2 == 0) {printf("该数是偶数。
\n");} else {printf("该数是奇数。
\n");}}```4. 求输入数字的平均值这个程序用于求输入数字的平均值。
```c#include <stdio.h>int main() {int count, i;double num, sum = 0.0, average;printf("请输入数字的个数:");scanf("%d", &count);printf("请输入这 %d 个数字:\n", count); for (i = 0; i < count; i++) {scanf("%lf", &num);sum += num;}average = sum / count;printf("平均值为:%lf\n", average);}```5. 判断闰年这个程序用于判断一个年份是否为闰年。
鲍威尔法实验报告姓名:刘阳阳学号:1201013009鲍威尔法程序设计1 方法(算法)介绍1.1 方法原理简介鲍威尔法——多维无约束优化算法是在无约束优化算法之一,首先选取一组共轭方向,从某个初始点出发,求目标函数在这些方向上的极小值点,然后以该点为新的出发点,重复这一过程直到获得满意解,其优点是不必计算目标函数的梯度就可以在有限步内找到极值点。
鲍威尔法是以共轭方向为基础的收敛较快的直接法之一,是一种十分有效的算法。
在无约束方法中许多算法都是以共轭方向作为搜索方向,它们具有许多特点。
根据构造共轭方向的原理不同,可以形成不同的共轭方向法。
基本算法1)任选一初始点x0,再选两个线性无关的向量,如坐标轴单位向量e1=[1,0]T 和e2=[0,1]T作为初始搜索方向。
2)从x0出发,顺次沿e1,e2作一维搜索,得x10,x20点,两点连线得一新方向d1=x20-x0。
用d1代替e1形成两个线性无关向量d1 ,e2 ,作为下一轮迭代的搜索方向。
再x20 出发,沿d1 作一维搜索得点x01,作为下一轮迭代的初始点。
3)从x1 出发,顺次沿,e2。
d1 作一维搜索,得到点x11,x21,两点连线得一新方向:d2=x21-x11。
4)沿d2d2作一维搜索得点.x2 ,即是二维问题的极小点x*把二维情况的基本算法扩展到n维,则鲍威尔基本算法的要点是:在每一轮迭代中总有一个始点(第一轮的始点是任选的初始点)和n个线性独立的搜索方向。
从始点出发顺次沿n个方向作一维搜索得一终点,由始点和终点决定了一个新的搜索方向。
用这个方向替换原来n个方向中的一个,于是形成新的搜索方向组。
替换的原则是去掉原方向组的第一个方向而将新方向排在原方向的最后。
此外规定,从这一轮的搜索终点出发沿新的搜索方向作一维搜索而得到的极小点,作为下一轮迭代的始点。
这样就形成算法的循环。
1.2 程序设计框图1.3程序说明本程序使用VB程序语言设计,在该程序中插入command模块picture模块利用矩阵函数语句迭代计算出鲍威尔最小值并将之在picture模块中显示出2 例题计算2.1 例题:用鲍威尔法求f(x1,x2)=10*(x1+x2-5)^2+(x1-x2)^22.2 计算结果与分析VB程序计算结果如图:经演算迭代后得到结果为:X0^2= 2.49952.5091 F0=f0(x0^2)=0.0008可见已足够接近极值点x*=(2.5,2.5)T 及极小值f(x*)=02.3 结论经验算鲍威尔法VB程序计算结果正确且为最优解。
#include <iostream.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <malloc.h>
#include <conio.h>
#define n1 2
#define tt 0.005
#define ad 0.0000001
//定义常量
//tt初始迭代步长
//ad收敛精度
float ia;
float fny(float *x)
{
float f;
f=10*pow((x[0]+x[1]-5),2)+pow((x[0]-x[1]),2); //目标函数return(f);
}
float *iterate(float *x,float a, float *s)
{
float *x1;
x1=(float *)malloc(n1 * sizeof(float));
for (int i=0;i<n1;i++)
x1[i]=x[i]+a*s[i];
return(x1);
}
float func(float *x,float a,float *s)
{
float *x1;
x1=iterate(x,a,s);
float f=fny(x1);
return(f);
}
void finding(float a[3],float f[3],float *xk,float *s) {
float t=tt;
float a1,f1;
a[0]=0; f[0]=func(xk,a[0],s);
for (int i=0;;i++)
{
a[1]=a[0]+t; f[1]=func(xk,a[1],s);
if (f[1]<f[0])
break;
if (fabs(f[1]-f[0])>=ad)
{
t=-t;
a[0]=a[1]; f[0]=f[1];
}
else{
if (ia==1) return;
t=t/2; ia=1;
}
}
for (i=0;;i++)
{
a[2]=a[1]+t; f[2]=func(xk,a[2],s);
if (f[2]>f[1])
break;
t=2*t;
a[0]=a[1]; f[0]=f[1];
a[1]=a[2]; f[1]=f[2];
}
if (a[0]>a[2])
{
a1=a[0]; f1=f[0];
a[0]=a[2]; f[0]=f[2];
a[2]=a1; f[2]=f1;
}
return;
}
//second insert
float lagrange(float *xk,float *ft,float *s)
{
float a[3],f[3];
float b,c,d,aa;
finding (a,f,xk,s);
for (int i=0;;i++)
{
if (ia==1)
{
aa=a[1]; *ft=f[1];
break;
}
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);
}
}
}
if (*ft>f[1])
{
aa=a[1]; *ft=f[1];
}
return (aa);
}
float *powell(float *xk)
{
float h[n1][n1],s[n1]={0,0},ff[n1+1]={0,0,0};
float f1,f3,aa;
float dk[n1],*x0,xk1[n1];
int m=0,i,j;
for (i=0;i<n1;i++)
{
for(j=0;j<n1;j++)
{
h[i][j]=0;
if (j==i)
h[i][j]=1;
}
}
for (int k=0;;k++)
{
ff[0]=fny(xk);
x0=xk;
for (i=0;i<n1;i++)
{
for (j=0;j<n1;j++)
s[j]=h[i][j];
float aa=lagrange(xk,&ff[i+1],s);
xk=iterate(xk,aa,s);
}
for (i=0;i<n1;i++)
{
float a,b;
dk[i]=xk[i]-x0[i];
xk1[i]=2*xk[i]-x0[i];
}
float max=fabs(ff[1]-ff[0]);
for (i=1;i<n1;i++)
if (fabs(ff[i+1]-ff[i])>max)
{
max=fabs(ff[i+1]-ff[i]);
m=i;
}
f3=fny(xk1);
if
((f3<ff[0])&&((ff[0]+f3-2*ff[2])*pow((ff[0]-ff[n1]-max),2)<0.5*max*pow((ff[0]-f3),2))) {
aa=lagrange(xk,&f1,dk);
xk=iterate(xk,aa,dk);
for (i=m;i<n1-1;i++)
for (j=0;j<n1;j++)
h[i][j]=h[i+1][j];
for (j=0;j<n1;j++)
h[n1-1][j]=dk[j];
}
else
{
if(ff[n1]>=f3) xk=xk1;
}
float xq=0;
for (i=0;i<n1;i++)
xq+=pow((xk[i]-x0[i]),2);
if (xq<=ad)
break;
}
return(xk);
}
void main()
{
float xk[n1]={0,0};//取初始点
float *xx;
xx=(float *)malloc(n1 *sizeof(float));
xx=powell(xk);
float ff=fny(xx);
cout<<"优化的结果为:"<<endl;
printf("\n\nThe Optimal Design Result Is:\n");
for (int i=0;i<n1;i++)
printf("\n\t x[%d] *=%f",i+1,xx[i]);
printf("\n\t f*=%f",ff);
getch();
}。