五.鮑威爾法－－已經調試好！！！
程序名PO
功能   解n維（n>=2)無約束優化問題
變量含義：
H0－－搜索區間步長
E1－－黃金分割精度
EP－－鮑威爾精度
N－－維數
X01  X02設計變量
有“!!!!”處語句為用戶應該修改的部分

主程序
80 REM MAIN PROGRAM:MPOW
81 PRINT"INPUT:N,E1,EP,H0"
82 INPUT N,E1,EP,H0
83 PRINT "INPUT:X01,X02" :REM !!!!!!!!!!!!!!!!!!!
84 DIM X(N),X0(N),XX1(N),XX2(N),XX3(N),S(N),SS(N,N+1),A(N),B(N),XP1(N),XP2(N),XK(N)
85 INPUT X0(1),X0(2):REM !!!!!!!!!!!!!!!!!!!!!!!!!
86 PRINT "PRIMAL DATA"
88 PRINT "N=";N,"E1=",E1,"EP=";EP,"H0=";H0,"X01=";X0(1),"X02=";X0(2):REM !!!!!!!!!!!!!!!
89 K=0:REM 循環次數
90 GOSUB 1000
92 PRINT"***************************"
94 PRINT "OPTIMUM RESULTS"
96 PRINT "X1=";X(1),"X2=";X(2),"F=";F,"K=";K:REM !!!!!!!!!!!!!!!!!!!!!!
98 END 




1000 REM FILE NAME:PO
1002 REM POWELLS METHOD
1006 FOR I=1 TO N
1008 FOR J=1 TO N+1
1010 SS(I,J)=0
1012 NEXT J
1014 SS(I,I)=1
1016 NEXT I
1018 FOR I=1 TO N
1020 XX1(I)=X0(I):X(I)=X0(I)
1022 NEXT I:GOSUB 8000:F1=FX:F0=F1
1024 DLT=-1
1026 FOR J=1 TO N
1028 FOR I=1 TO N
1030 S(I)=SS(I,J):X0(I)=X(I)
1031 NEXT I:GOSUB 1040
1032 DF=F0-F:F0=F
1033 IF DF<=DLT GOTO 1035
1034 DLT=DF:M=J
1035 NEXT J:GOTO 2400


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




2400 SDX=0
2402 FOR I=1 TO N
2404 SDX=SDX+(X(I)-XX1(I))^2
2406 NEXT I

2407 K=K+1:REM 循環次數
2408 IF SDX<=EP GOTO 2510
2410 FOR I=1 TO N
2412 XX2(I)=X(I)
2414 NEXT I:F2=F



2422 FOR I=1 TO N
2444 SS(I,N+1)=X(I)-XX1(I)
2446 S(I)=SS(I,N+1)
2447 X0(I)=X(I)
2448 NEXT I:GOSUB 1040
2466 FOR I=1 TO N
2467 XK(I)=X(I)
2468 XX3(I)=2*XX2(I)-XX1(I):X(I)=XX3(I)
2470 NEXT I:GOSUB 8000:F3=FX
2472 Q=(F1-2*F2+F3)*(F1-F2-DLT)^2
2474 D=0.5*DLT*(F1-F3)^2
2476 IF F3>=F1 GOTO 2496
2478 IF Q>=D GOTO 2496
2480 FOR J=M+1 TO N+1
2482 FOR I=1 TO N
2484 SS(I,J-1)=SS(I,J)
2486 NEXT I
2488 NEXT J
2490 FOR I=1 TO N
2492 X0(I)=XK(I)
2494 NEXT I:GOTO 1018
2496 IF F2>=F3 GOTO 2504 
2498 FOR I=1 TO N
2500 X0(I)=XX2(I)
2502 NEXT I:GOTO 1018
2504 FOR I=1 TO N
2506 X0(I)=XX3(I)
2508 NEXT I :GOTO 1018
2510 RETURN



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