[分享]dfp变度量法优化程序
<P>PROGRAM TON<BR>DIMENSION X(2),Y(2),GX(2)<BR>IOBF=1<BR>IDER=1<BR>X(1)=-1.2<BR>X(2)=1.5<BR>CALL VARMT1(X,Y,FY,GX,2,5000,ICOVG,0.000001,1)<BR>WRITE(6,1)Y,FY,ICOVG<BR>1 FORMAT(5X,3E21.10,2X,I5//)<BR>STOP<BR>END<BR><BR>SUBROUTINE VARMT1(X,Y,FY,GX,N,IMOS,ICOVG,TL,IT)<BR>DIMENSION X(N),Y(N),GX(N),H(20,20),DELX(20),DELG(20)<BR>IF(IT.GT.0) GO TO 151<BR>TL=0.0001<BR>151 ITER=0<BR>ICOVG=0<BR>ER1=0.01<BR>ER2=0.001<BR>ET1=0.001<BR>ET2=10.<BR>IRESET=N+1<BR>INDEX=IRESET<BR>CALL OBF(X,FX)<BR>CALL DER(X,GX)<BR>GXN=0.0<BR>DO 2 I=1,N<BR>DELG(I)=GX(I)<BR>2 GXN=GXN+GX(I)*GX(I)<BR>GXN=SQRT(GXN)<BR>1 IF(GXN.LE.TL) GO TO 60<BR>IF(ITER.GE.IMOS) GO TO 80<BR>IF(INDEX.NE.IRESET) GO TO 90<BR>5 INDEX=0<BR>DO 20 I=1,N<BR>DO 10 J=1,N<BR>10 H(I,J)=0.0<BR>H(I,I)=1.0<BR>20 GX(I)=-DELG(I)<BR>130 TD=ET2*GXN<BR>IF(TD.LT.ET1) GO TO 21<BR>TD=ET1<BR>21 CALL DSCPOW(X,FX,Y,FY,GX,GX,N,2,TD,1,ITER,IEX)<BR>IF(IEX.NE.1) GO TO 31<BR>ICOVG=-1<BR>RETURN<BR>31 ITER=ITER+1<BR>INDEX=INDEX+1<BR>CALL DER(Y,GX)<BR>GXN=0.0<BR>DO 50 I=1,N<BR>GXN=GXN+GX(I)*GX(I)<BR>DELG(I)=GX(I)-DELG(I)<BR>DELX(I)=Y(I)-X(I)<BR>IF(INDEX.NE.IRESET) GO TO 50<BR>DELG(I)=GX(I)<BR>50 X(I)=Y(I)<BR>FX=FY<BR>GXN=SQRT(GXN)<BR>GO TO 1<BR>90 CALL HDFP(H,N,DELG,DELX,IE)<BR>DO 120 I=1,N<BR>120 DELG(I)=GX(I)<BR>IF(IE.EQ.-1) GO TO 5<BR>SN=0.0<BR>GS=0.0<BR>DO 100 I=1,N<BR>GX(I)=0.0<BR>DO 110 J=1,N<BR>110 GX(I)=GX(I)-H(I,J)*DELG(J)<BR>GS=GS+GX(I)*DELG(I)<BR>100 SN=SN+GX(I)*GX(I)<BR>SN=SQRT(SN)<BR>ER=ER2*GXN<BR>IF(ER1.GT.ER) GO TO 52<BR>ER=ER1<BR>52 IF(GS.GT.-ER*GXN*SN) GO TO 5<BR>GO TO 130<BR>60 ICOVG=1<BR>80 RETURN<BR>END</P><P>SUBROUTINE DSCPOW(X,FX,Y,FY,S,DELX,N,INDIC,TOL,ITOL,ITER,IEX)<BR>DIMENSION X(N),Y(N),S(N),DELX(N)</P>
<P>IF(ITOL.EQ.1)GO TO 110<BR>TOL=.001<BR>110 FTOL2=TOL/100<BR>IEX=1<BR>K=-2<BR>M=0<BR>FA=FX<BR>DA=0.<BR>STEP=1.0<BR>D=STEP<BR>IF(INDIC.EQ.2)GO TO 1<BR>IF(ITER.EQ.0)GO TO 1<BR>DXNORM=0.<BR>SNORM=0.<BR>DO 102 I=1,N<BR>DXNORM=DXNORM+DELX(I)*DELX(I)<BR>102 SNORM=SNORM+S(I)*S(I)<BR>IF(INDIC.NE.1)GO TO 103<BR>IF(DXNORM.GE.SNORM)GO TO 1<BR>103 RATIO=DXNORM/SNORM<BR>STEP=SQRT(RATIO)<BR>D=STEP<BR>1 DO 2 I=1,N<BR>2 Y(I)=X(I)+D*S(I)<BR>CALL OBF(Y,F)<BR>K=K+1<BR>IF(F-FA)5,3,6<BR>3 D=0.5*(DA+D)<BR>DO 4 I=1,N<BR>4 Y(I)=X(I)+D*S(I)<BR>CALL OBF(Y,F)<BR>IF(F-FA)204,202,205<BR>202 DO 203 I=1,N<BR>203 Y(I)=X(I)+DA*S(I)<BR>FY=FA<BR>GO TO 326<BR>204 FC=FA<BR>DC=DA+2.*(D-DA)<BR>FB=F<BR>DB=D<BR>GO TO 21<BR>205 STEP=0.5*STEP<BR>IF(K)7,8,206<BR>206 DC=D<BR>D=DA<BR>DA=DB<BR>DB=D<BR>FC=F<BR>F=FA<BR>FA=FB<BR>FB=F<BR>GO TO 10<BR>5 FC=FB<BR>FB=FA<BR>FA=F<BR>DC=DB<BR>DB=DA<BR>DA=D<BR>D=2.0*D+STEP<BR>GO TO 1<BR>6 IF(k)7,8,9<BR>7 FB=F<BR>DB=D<BR>D=-D<BR>STEP=-STEP<BR>GO TO 1<BR>8 FC=FB<BR>FB=FA<BR>FA=F<BR>DC=DB<BR>DB=DA<BR>DA=D<BR>GO TO 21<BR>9 DC=DB<BR>DB=DA<BR>DA=D<BR>FC=FB<BR>FB=FA<BR>FA=F<BR>10 D=0.5*(DA+DB)<BR>DO 11 I=1,N<BR>11 Y(I)=X(I)+D*S(I)<BR>CALL OBF(Y,F)<BR>12 IF((DC-D)*(D-DB))15,13,18<BR>13 DO 14 I=1,N<BR>14 Y(I)=X(I)+DB*S(I)<BR>FY=FB<BR>IF(IEX.EQ.0)GO TO 32<BR>GO TO 326<BR>15 IF(F-FB)16,13,17<BR>16 FC=FB<BR>FB=F<BR>DC=DB<BR>DB=D<BR>GO TO 21<BR>17 FA=F<BR>DA=D<BR>GO TO 21<BR>18 IF(F-FB)19,13,20<BR>19 FA=FB<BR>FB=F<BR>DA=DB<BR>DB=D<BR>GO TO 21<BR>20 FC=F<BR>DC=D<BR>21 IF(ABS(DA-DB).GE.ABS(DC-DB))GO TO 50<BR>D=DA<BR>DA=DC<BR>DC=D<BR>F=FA<BR>FA=FC<BR>FC=F<BR>50 IF(ABS((DC-DB)/(DA-DB)).GT.0.25)GO TO 26<BR>D=0.5*(DA+DB)<BR>DO 51 I=1,N<BR>51 Y(I)=X(I)+D*S(I)<BR>CALL OBF(Y,F)<BR>IF(F-FB)52,13,53<BR>52 FC=FB<BR>FB=F<BR>DC=DB<BR>DB=D<BR>GO TO 26<BR>53 FA=F<BR>DA=D<BR>26 A=FA*(DB-DC)+FB*(DC-DA)+FC*(DA-DB)<BR>IF(A)22,30,22<BR>22 D=0.5*((DB*DB-DC*DC)*FA+(DC*DC-DA*DA)*FB+(DA*DA-DB*DB)*FC)/A<BR>IF((DA-D)*(D-DC))13,13,23<BR>23 DO 24 I=1,N<BR>24 Y(I)=X(I)+D*S(I)<BR>CALL OBF(Y,F)<BR>IF(ABS(FB-F)-(ABS(FB)+FTOL2)*TOL)28,28,12<BR>28 IEX=0<BR>IF(F-FB)29,13,13<BR>29 FY=F<BR>GO TO 32<BR>30 IF(M)31,31,13<BR>31 M=M+1<BR>GO TO 10<BR>32 DO 99 I=1,N<BR>IF(Y(I).NE.X(I))GO TO 326<BR>99 CONTINUE<BR>33 IF(NTOL.EQ.5)GO TO 34<BR>IEX=1<BR>NTOL=NTOL+1<BR>TOL=TOL/10.<BR>GO TO 12<BR>34 IEX=1<BR>GO TO 35<BR>326 IF(FY.LT.FX)GO TO 36<BR>IEX=1<BR>GO TO 35<BR>36 IF(IEX.EQ.0)GO TO 35<BR>IEX=2<BR>35 RETURN<BR>END</P>
<P><BR><BR>SUBROUTINE OBF(X,FX)<BR>DIMENSION X(2)<BR>FX=100*(X(2)-X(1)*X(1))*(X(2)-X(1)*X(1))+(1-X(1))*(1-X(1))<BR>RETURN<BR>END</P>
<P>SUBROUTINE HDFP(H,N,DELG,DELX,IE)<BR>DIMENSION H(20,20),DELX(20),DELG(20),DHG(20),DGH(20)<BR>DXDG=0.0<BR>DGHDG=0.0<BR>DO 20 I=1,N<BR>DHG(I)=0.0<BR>DGH(I)=0.0<BR>DO 10 J=1,N<BR>DHG(I)=DHG(I)-H(I,J)*DELG(J)<BR>10 DGH(I)=DGH(I)+DELG(J)*H(J,I)<BR>DXDG=DXDG+DELX(I)*DELG(I)<BR>20 DGHDG=DGHDG+DGH(I)*DELG(I)<BR>IE=-1<BR>IF(ABS(DXDG).LT.0.1E-60) GO TO 40<BR>IF(DGHDG.LT.0.1E-60) GO TO 40<BR>IE=0<BR>DO 30 I=1,N<BR>DO 30 J=1,N<BR>30 H(I,J)=H(I,J)+DELX(I)*DELX(J)/DXDG+DHG(I)*DGH(J)/DGHDG<BR>40 RETURN<BR>END<BR><BR><BR>SUBROUTINE DER(X,GX)<BR>DIMENSION X(2),GX(2)<BR>GX(1)=-400*X(1)*(X(2)-X(1)*X(1))-2*(1-X(1))<BR>GX(2)=200*(X(2)-X(1)*X(1))<BR>RETURN<BR>END</P> fortran 编的啊，看不懂，楼主有没有VB的dfp程序代码 <P>谢谢了，请问一下楼主，这是不是就是变尺度法</P> <P>请问一下楼主， 我要用fortran编一个用DFP法算的电力系统潮流计算程序，你能把这个程序说明一下吗，我好借用一下？？我的邮箱是<a href="mailtxuang34@163.com" target="_blank" >xuang34@163.com</A>,谢谢了!!</P> 谢谢楼主,请问有没有BFGS程序呀 <P>我发现这个网站真是太好!!<BR></P>
页:
[1]