罚函数法MATLAB程序

一、进退法、0.618法、Powell法、罚函数法的Matlab程序设计罚函数法（通用）function y=ff(x,k)y=-17.86*0.42*x(1)/(0.8+0.42*x(1)*(1-exp(-2*(0.8+0.42*x(1)/3)*exp(-1.6)*x(2)-22.99*x(1)/(0.8+x(1)*(1-exp(-2*(0.8+x(1)/3)*x(3)+k*(x(2)-(1.22*102*(9517.8*exp(-1.6-2*0.42*x(1)/3)*x(2)+19035.6*exp(-2*x(1)/3)*x(3)/(1.22*102+9517.8*exp(-1.6-2*0.42*x(1)/3)*x(2)+19035.6*exp(-2*x(1)/3)*x(3)2+k*(x(3)-exp(-0.8-2*x(1)/3)*x(3)-exp(-2.4-2*0.42*x(1)/3)*x(2)2;% 主函数，参数包括未知数的个数n，惩罚因子q，惩罚因子增长系数k，初值x0，以及允许的误差rfunction G=FHS(x0,q,k,n,r,h,a) l=1;while (l) x=powell(x0,n,q,r(1),h,a); %调用powell函数 g(1)=ff1(x),g(2)=ff2(x) . . . g(p)=ffp(x); %调用不等式约束函数，将其值 %存入数组g h(1)=hh1(x),h(2)=hh2(x) . . . h(t)=hht(x); %调用等式约束函数，将其值%存入数组h for i=1:p if g(i)r for j=1:t if abs(h(j)r(2) continue; else break; end end else break; end end if (i=p)&(j=t) %如果所有约束条件的值都在允许的范围内就结束%搜索，否则加大惩罚力度继续搜索 l=0; else x0=x; q=k*q; endend G=x%powell算法，用于寻找无约束最优值点 function powel=powell(x0,n,q,r,h,a) d=eye(n); %n个线性无关的初始搜索方向 k=1; kk=1; xx(1,1:n)=x0; while (kk) y(1,1:n)=xx(k,1:n); for j=1:n s(j)=HJ(y(j,1:n),d(j,1:n),q); %调用0.618算法 y(j+1,1:n)=y(j,1:n)+s(j).*d(j,1:n); end d(n+1,1:n)=y(n+1,1:n)-y(1,1:n); if (norm(d(n+1,1:n),2)=ww) m=i; end end cha=ff(y(1,1:n),q)-2*ff(y(n+1,1:n),q)+ff(2*y(n+1,1:n)-y(1,1:n),q); cha1=2*(ff(y(m,1:n),q)-ff(y(m+1,1:n),q); if (cha=0.1) if y1y2 b=r2; r2=r1; y2=y1; r1=a+(1-r)*(b-a); y1=ff(x0+r1.*dd,q); else a=r1; r1=r2; y1=y2; r2=a+r*(b-a); y2=ff(x0+r2.*dd,q); end end w=(r1+r2)/2%进退法function a,b=JTF(x0,h,d,dd,q) r0=0; y0=ff(x0+r0.*dd,q); k=0; l=1; while (l) r1=r0+h; y1=ff(x0+r1.*dd,q); if y1y0 h=d*h; r=r0; r0=r1; y0=y1; else if k=0; h=-h; r=r0; else l=0; break; end end k=k+1; end a=min(r,r1); b=max(r,r1);二、进退法、0.618法的C语言程序设计例1 设f(x)=x2-2*x+4 ，试确定初始搜索区间，并输出迭代次数main（）int n=0；float t，h，a，b，c，d，f1，f2；printf（“t，h：”）；scanf（“%f，%f”，&t，&h）；f1=t*t-2*t+4；if（t+h）*（t+h）-2*（t+h）+4=f1）c=t+h;f2=c*c-2*c+4;elseh=-h;c=t+h; f2=c*c-2*c+4;while(f2=f1)h=2*h;d=t;t=c;f1=f2;c=t+h; f2=c*c-2*c+4;n=n+1;if(c=d)a=c;b=d;elsea=d;b=c;printf(“a,b=%.2f,%fn”,a,b);printf(“n=%dn”,n);运行情况如下：（1）t,h:0,0.05 a,b=0.35,1.55 n=4(2) t,h:-1.0,0.02 a,b=0.26,4.10 n=7(3)t,h:-100,0.01 a,b=-59.05,63.83 n=13例2 用0.618法求f(x)=8x3-2*x2-7*x+3 的局部最优解允许误差=0.0001 ，初始点设为x0=0 ，初始步长h=1 #include“math.h”float f(x)float x；float y;y=(8.0*x-2.0)*x-7.0)*x+3.0;return(y);main()int n=0;float t,h,a,b,c,d,x0,f0,f1,f2,fc,fd;printf(“t,h:”);scanf(“%f,%f”,&t,&h);f1=f(t);if (f(t+h)=f1)c=t+h;elseh=-h;c=t+h;f2=f(c);while(f2=f1)d=t;t=c;f1=f2;c=t+h;f2=f(c);if(c=fd) a=c;c=d;fc=fd;d=a+0.618*(b-a); fd=f(d);else b=d;d=c;fd=fc;c=a+0.382*(b-a);f