七.DFP變尺度法
程序名 PDF
功能 解無約束優化
主程序
######################
輸入變量含義：
H0－－搜索區間步長
E1－－黃金分割精度
ED－－DFP精度
N－－維數
X01  X02設計變量
有“!!!!”處為用戶應該修改的部分




300 REM MAIN PROGRAM OF DFP METHOD
301 PRINT "INPUT:N,H0,E1,ED"
302 INPUT N,H0,E1,ED
304 DIM X(N),X0(N),SI(N),Y(N),T(N),S(N),R(N),AA(N,N),C(N,N),D(N,N),HH(N,N),U(N,N),M(N,N),P(N,N)
305 PRINT "INPUT:X01,X02":REM !!!!!!!!!!!!!!!!!!!!!!!!!!
306 INPUT X0(1),X0(2):REM !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
310 PRINT "PRIMAL DATA"
312 PRINT "N=";N,"H0=";H0,"E1=";E1,"ED=";ED,"X01=";X0(1),"X02=";X0(2):REM !!
313 GOSUB 3000
314 PRINT "**********************************"
316 PRINT "OPTIMUM RESULTS"
318 PRINT "X1=";X(1),"X2=";X(2),"F=";F:REM !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
320 END



1000 REM FILE NAME :DFP
1020 REM SUBPROGRAM OF DFP

1040 REM DETERMINATION OF THE INITIAL INTERVAL黃金分割法
1041 H=H0
1042 FOR I=1 TO N
1044 X1(I)=X0(I):X(I)=X1(I)
1046 NEXT I
1048 GOSUB 8000:Y1=FX
1050 FOR I=1 TO N
1052 X2(I)=X1(I)+H*S(I):X(I)=X2(I)
1054 NEXT I
1056 GOSUB 8000:Y2=FX
1058 IF Y2<Y1 GOTO 1074
1060 H=-H0
1062 FOR I=1 TO N
1064 X3(I)=X1(I)
1066 NEXT I:Y3=Y1
1068 FOR I=1 TO N
1070 X1(I)=X2(I):X2(I)=X3(I)
1072 NEXT I:Y1=Y2:Y2=Y3
1074 H=2*H
1076 FOR I=1 TO N
1078 X3(I)=X2(I)+H*S(I):X(I)=X3(I)
1080 NEXT I:GOSUB 8000:Y3=FX
1082 IF Y2>=Y3 GOTO 1068
1084 IF H<0 GOTO 1092
1086 FOR I=1 TO N
1088 A(I)=X1(I):B(I)=X3(I)
1090 NEXT I:GOTO 1400
1092 FOR I=1 TO N
1094 A(I)=X3(I):B(I)=X1(I)
1096 NEXT I
1400 REM FILE NAME:GS
1402 REM GOLDEN SECTION METHOD
1406 FOR I=1 TO N
1408 XP1(I)=A(I)+0.382*(B(I)-A(I))
1409 X(I)=XP1(I)
1410 NEXT I:GOSUB 8000:Y1=FX
1412 FOR I=1 TO N
1414 XP2(I)=A(I)+0.618*(B(I)-A(I))
1415 X(I)=XP2(I) 
1416 NEXT I:GOSUB 8000:Y2=FX
1418 IF Y1<Y2 GOTO 1432
1420 FOR I=1 TO N
1422 A(I)=XP1(I):XP1(I)=XP2(I)
1424 NEXT I:Y1=Y2
1426 FOR I=1 TO N
1428 XP2(I)=A(I)+0.618*(B(I)-A(I))
1429 X(I)=XP2(I)
1430 NEXT I:GOSUB 8000:Y2=FX:GOTO 1444
1432 FOR I=1 TO N
1434 B(I)=XP2(I):XP2(I)=XP1(I)
1436 NEXT I:Y2=Y1
1438 FOR I=1 TO N
1440 XP1(I)=A(I)+0.382*(B(I)-A(I))
1441 X(I)=XP1(I):NEXT I
1442 GOSUB 8000:Y1=FX
1444 L=0
1446 FOR I=1 TO N
1448 L=L+(B(I)-A(I))^2
1450 NEXT I
1452 L=SQR(L)
1454 IF L>E1 GOTO 1418 
1456 FOR I=1 TO N
1458 X(I)=0.5*(B(I)+A(I))
1460 NEXT I :GOSUB 8000:F=FX
1462 RETURN：REM 黃金分割法
 

3000 REM FILE NAME :DFP
3002 REM DFP METHOD
3008 FOR I=1 TO N 
3010 X(I)=X0(I)
3011 NEXT I
3020 GOSUB 8100:GOSUB 8000
3022 K=1
3024 FOR I=1 TO N
3026 FOR J=1 TO N
3028 AA(I,J)=0
3030 NEXT J
3032 AA(I,I)=1
3034 NEXT I
3035 FOR I=1 TO N
3036 S(I)=-T(I)
3037 NEXT I
3038 Q=0
3039 FOR I=1 TO N
2040 Q=Q+T(I)^2
3041 NEXT I
3042 Q=SQR(Q)
3048 IF Q<=ED GOTO 3170 
3049 IF K=1 GOTO 3140
3050 FOR I=1 TO N
3052 SI(I)=X(I)-SI(I):Y(I)=T(I)-Y(I)
3054 NEXT I
3056 FOR I=1 TO N
3058 FOR J=1 TO N
3060 C(I,J)=SI(I)*SI(J)
3062 D(I,J)=Y(I)*Y(J)
3064 NEXT J
3066 NEXT I
3068 FOR I=1 TO N
3070 FOR J=1 TO N
3072 HH(I,J)=0
3074 FOR L=1 TO N
3076 HH(I,J)=HH(I,J)+AA(I,L)*D(L,J)
3078 NEXT L
3080 NEXT J
3082 NEXT I
3084 FOR I= 1 TO N
3086 FOR J=1 TO N
3088 U(I,J)=0
3090 FOR L=1 TO N
3092 U(I,J)=U(I,J)+HH(I,L)*AA(L,J)
3093 NEXT L:NEXT J:NEXT I
3094 W1=0:W2=0
3096 FOR I=1 TO N
3098 W1=W1+SI(I)*Y(I)
3099 NEXT I
3100 FOR I=1 TO N
3101 R(I)=0
3102 FOR J=1 TO N
3106 R(I)=R(I)+Y(J)*AA(J,I)
3108 NEXT J
3110 NEXT I
3112 FOR I=1 TO N
3114 W2=W2+R(I)*Y(I)
3116 NEXT I
3118 FOR I=1 TO N
3120 FOR J=1 TO N
3122 M(I,J)=C(I,J)/W1
3124 P(I,J)=U(I,J)/W2
3126 AA(I,J)=AA(I,J)+M(I,J)-P(I,J)
3128 NEXT J
3130 NEXT I
3132 FOR I=1 TO N
3133 S(I)=0
3134 FOR J=1 TO N
3136 S(I)=S(I)-AA(I,J)*T(J)
3138 NEXT J:NEXT I
3140 FOR I=1 TO N
3142 SI(I)=X(I):Y(I)=T(I)
3143 NEXT I
3144 FOR I=1 TO N
3145 X0(I)=X(I)
3146 NEXT I
3147 GOSUB 1040:GOSUB 8100
3150 K=K+1 
3152 IF K>N GOTO 3156
3154 GOTO 3038
3156 Q=0
3158 FOR I=1 TO N
3160 Q=Q+T(I)^2
3162 NEXT I
3163 Q=SQR(Q)
3164 IF Q<=ED GOTO 3170
3166 GOTO 3022
3170 RETURN



8000 REM OBJECTIVE FUNCTION
8001 FX=X(1)^2+X(2)^2-X(1)*X(2)-10*X(1)-4*X(2)+60
8002 REM !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
8004 RETURN

8100 REM FIVE POINT CENTRAL FORMAL
8102 FOR I=1 TO N
8104 IF X(I)=0 GOTO 8108
8106 TT=X(I)/100:GOTO 8110
8108 TT=0.01
8110 FOR J=1 TO 5
8112 IF J=3 GOTO 8120
8114 X(I)=X(I)+(J-3)*TT
8116 GOSUB 8000:Z(J)=FX
8118 X(I)=X(I)-(J-3)*TT
8120 NEXT J
8122 T(I)=(Z(1)+8*Z(4)-8*Z(2)-Z(5))/12/TT
8124 NEXT I
8126 RETURN