罰函數(shù)法MATLAB程序

1、一、進(jìn)退法、0.618法、Powell法、罰函數(shù)法的 Matlab程序設(shè)計罰函數(shù)法(通用)fun ctio n 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*10A2*(9517.8*exp(-1.6-2*0.42*x(1)/3)*x(2)+19035.6*exp(-2*x(1)/3)*x(3)/(1.22*10A2+9517.8*exp(-

2、1.6-2*0.42*x(1)/3)*x(2)+19035.6*exp(-2*x(1)/3)*x(3)A2+k*(x(3)-exp(-0.8-2*x(1)/3)*x(3) -exp(-2.4-2*0.42*x(1)/3)*x(2)A2;%主函數(shù)，參數(shù)包括未知數(shù)的個數(shù)n,懲罰因子q，懲罰因子增長系數(shù) k,初值x0,以及允許的誤差rfun ctio n G=FHS(x0,q,k, n,r,h,a)l=1;while (l)x=powell(x0,n,q,r(1),h,a);% 調(diào)用 powell 函數(shù)g(1)=ff1(x),g(2)=ff2(x) . . . g(p)=ffp(x);% 調(diào)用不等式

3、約束函數(shù)，將其值%存入數(shù)組gh(1)=hh1(x),h(2)=hh2(x) . . . h(t)=hht(x);% 調(diào)用等式約束函數(shù)，將其值%存入數(shù)組hfor i=1:pif g(i)<rfor j=1:tif abs(h(j)<r(2)con ti nue;elsebreak;endendelsebreak;endendif (i=p)&(j=t)%如果所有約束條件的值都在允許的范圍內(nèi)就結(jié)束%搜索，否則加大懲罰力度繼續(xù)搜索l=0;elsex0=x;q=k*q;endendG=x%powell算法，用于尋找無約束最優(yōu)值點fun cti onpowel=powell(x0 ,

4、n, q,r,h,a)d=eye(n);%n個線性無關(guān)的初始搜索方向k=1;kk=1;xx(1,1: n)=xO;while (kk)y(1,1: n)=xx(k,1: n);for j=1: ns(j)=HJ(y(j,1:n),d(j,1:n),q);% 調(diào)用 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 (n orm(d( n+1,1: n),2)<r)kk=0;break;elseww=0;m=1;for i=1: ngg=ff(y(i,1: n),q)

5、-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); cha仁 2*(ff(y(m,1: n),q)-ff(y(m+1,1: n),q);if (cha<cha1)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) for j=m+1: nd(j,1: n)=d(j+1,1: n);e

6、ndk=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);%調(diào)用進(jìn)退法算法，確定范圍r=0.618;r1= a+(1-r)*(b-a);r2=a+r*(b-a);y仁 ff(x0+r1.*dd,q);y2=ff(x0+r2.*dd,q);k=1;while (abs(r1-r2)>=0.1)if y1<y2b=r2; r2=r1;y2=y1; r1=a+(1-r)*(b-a); y1= ff

7、(xO+r1.*dd,q);elsea=r1;r1=r2; y1=y2; r2=a+r*(b-a); y2=ff(x0+r2.*dd,q);endendw=(r1+r2)/2%進(jìn)退法fun ctio n a,b=JTF(xO,h,d,dd,q) r0=0;yO=ff(xO+rO.*dd,q);k=0;l=1; while (l) r仁rO+h;y1= ff(x0+r1.*dd,q);if y1<y0 h=d*h; r=r0; r0=r1; y0=y1;elseif k=0;h=-h;r=r0;elsel=0;break;endendk=k+1;enda=min(r ,r1);b=max(

8、r ,r1);二、進(jìn)退法、0.618法的C語言程序設(shè)計例1設(shè)f(x)=x2-2*x+4，試確定初始搜索區(qū)間，并輸出迭代次數(shù)main ()int n=0 ;float t， h， a，b， c， d， fl， f2 ；printf ( “， 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

9、; 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.05a,b=0.35,1.55n=4(2) t,h:-1.0,0.02a,b=0.26,4.10n=7(3) t,h:-100,0.01a,b=-59.05,63.83n=13,初始點例2 用0.618法求f(x)=8xA3-2*xA2-7*x+3的局部最優(yōu)解.允許誤差8=0.0001設(shè)為x0=0 ，初始步長h=1.#include math.h

10、”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: ”)；seanf( %f,%f ",&t,&h);f仁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;f 1= f2;c=t+h;f2=f(c);if(c<=d)a=c;b=d;elsea=d;b=c; c=a+0.382*(b-a);fc=f(c);d=a+0.618*(b-a);fd=f(d); doif (fc>=fd) a=c;c=d;fc=fd;d=a+0.618*(b-a); fd=f(d);else b=d