罚函数法MATLAB程序

一、进退法、0.618法、Powell法、罚函数法的Matlab程序设计罚函数法（通用）functiony=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*10^2*(9517.8*exp(-1.6-2*0.42*x(1)/3)*x(2)+19035.6*exp(-2*x(1)/3)*x(3)))/(1.22*10^2+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，以及允许的误差rfunctionG=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);%调用不等式约束函数，将其值%存入数组gh(1)=hh1(x),h(2)=hh2(x)...h(t)=hht(x);%调用等式约束函数，将其值%存入数组hfori=1:pifg(i)rforj=1:tifabs(h(j))r(2)continue;elsebreak;endendelsebreak;endendif(i==p)&(j==t)%如果所有约束条件的值都在允许的范围内就结束%搜索，否则加大惩罚力度继续搜索l=0;elsex0=x;q=k*q;endendG=x%powell算法，用于寻找无约束最优值点functionpowel=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);forj=1:ns(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);endd(n+1,1:n)=y(n+1,1:n)-y(1,1:n);if(norm(d(n+1,1:n),2)r)kk=0;break;elseww=0;m=1;fori=1:ngg=ff(y(i,1:n),q)-ff(y(i+1,1:n),q);if(gg=ww)m=i;endendcha=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(chacha1)s(n+1)=HJ(y(n+1,1:n),h,a,d(n+1,1:n),q)xx(k+1,1:2)=y(n+1,1:n)+s(n+1).*d(n+1,1:n)forj=m+1:nd(j,1:n)=d(j+1,1:n);endk=k+1;elsexx(k+1,1:n)=y(n+1,1:n);k=k+1;endendendpowel=y(n+1,1:n)functionw=HJ(x0,h,d,dd,q)%0.618算法[a,b]=JTF(x0,h,d,dd,q);%调用进退法算法，确定范围r=0.618;r1=a+(1-r)*(b-a);r2=a+r*(b-a);y1=ff(x0+r1.*dd,q);y2=ff(x0+r2.*dd,q);k=1;while(abs(r1-r2)=0.1)ify1y2b=r2;r2=r1;y2=y1;r1=a+(1-r)*(b-a);y1=ff(x0+r1.*dd,q);elsea=r1;r1=r2;y1=y2;r2=a+r*(b-a);y2=ff(x0+r2.*dd,q);endendw=(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);ify1y0h=d*h;r=r0;r0=r1;y0=y1;elseifk==0;h=-h;r=r0;elsel=0;break;endendk=k+1;enda=min(r,r1);b=max(r,r1);二、进退法、0.618法的C语言程序设计例1设f(x)=x^2-2*x+4，试确定初始搜索区间，并输出迭代次数main（）{intn=0；floatt，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;}else{h=-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;}else{a=d;b=c;}printf(“[a,b]=[%.2f,%f]\n”,a,b);printf(“n=%d\n”,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)=8x^3-2*x^2-7*x+3的局部最优解．允许误差ε=0.0001，初始点设为x0=0，初始步长h=1．#include“math.h”floatf(x)floatx；floaty;{y=((8.0*x-2.0)*x-7.0)*x+3.0;return(y);}main(){intn=0;floatt,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;else{h=-h;c=t+h;}f2=f(c);while(f2=f1){d=t;t=c;f1=f2;c=t+h;f2=f(c);}if(c=d){a=c;b=d;}else{a=d;b=c;}c=a+0.382*(b-a);fc=f(c);d=a+0.618*(b-a);fd=f(d);do{if(fc=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);fc=f(c);}n++;}while(b-a0.0001);x0=(a+b)/2.0;f0=f(x0);pintf(“n=%d,x0=%7.4f,f0=%7.4f\n”,n,x0,f0);}运行情况如下：t,h:0,1.0n=21,xo=0.6298,f0=-0.2034．