工程优化中常用的几种方法和实现

调试通过
功能：求解n维具有不等式约束优化问题的最优解
变量含义：
N－－维数
KG－－约束数
H0－－搜索区间步长(0.1??)
R0－－初始惩罚因子(.001??)
C--惩罚系数（>1）
E1－－一维搜索法黄金分割精度(0.001??)
EP－－鲍威尔精度(0.001???)
ES--惩罚函数法精度
X01  X02设计变量
有“!!!!”处语句为用户应该修改的部分


子程序：


100 REM MAIN:MEP
101 PRINT "INPUT:N,KG,H0,R0,C,E1,EP,ES"
102 INPUT N,KG,H0,R0,C,E1,EP,ES
103 DIM X(N),X0(N),XX1(N),XX2(N),XX3(N),S(N),SS(N,N+1),A(N),B(N)
104 DIM XP1(N),XP2(N),G(KG),XK(N)
105 PRINT "INPUT:X01,X02"：INPUT X0(1),X0(2):REM !!!
106 PRINT "PRIMAL DATA"
107 PRINT "N=";N,"KG=";KG,"H0=";H0,"R0=";R0,"C=";C,"E1=";E1,"EP=";EP
108 PRINT "ES=";ES,"X01=";X0(1),"X02=";X0(2)
110 GOSUB 5200
114 PRINT "OPTIMUM RESULTS"
116 PRINT "X1=";X(1),"X2=";X(2),"F=";F:REM !!!
120 END


1000 REM FILE NAME:PO(1000--2510鲍威尔法）
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















5200 REM FILE NAME:EP
5206 R=R0
5208 K=1
5210 GOSUB 1000
5212 IF K=1 GOTO 5216 
5214 IF ABS((FU-F00)/FU)<=ES GOTO 
5216 K=K+1
5218 R=R*C
5220 FOR I=1 TO N
5222 X0(I)=X(I)
5224 NEXT I
5226 F00=FU
5228 GOTO 5210
5230 GOSUB 9000
5231 F=FU
5232 RETURN



8000 REM CALCULATION OF EXTERIOR PENALTY FUNTION VALUE
8002 FX=X(1)^2+X(2)^2-X(1)*X(2)-10*X(1)-4*X(2)+60
8003 GOSUB 9000
8004 SG=0
8005 FOR I=1 TO KG
8006 IF G(I)>=0 GOTO 8008
8007 GOTO 8012
8008 NEXT I
8009 GOTO 8017
8012 FOR I=1 TO KG
8014 SG=SG+G(I)^2
8016 NEXT I
8017 FX=FU+R*SG
8018 RETURN



9000 REM FUNCTION
9004 FU=(X(1)-8)^2+(X(2)-8)^2:REM !!!
9006 G(1)=X(1):REM !!!
9008 G(2)=X(2)-1:REM !!!
9014 G(3)=11-X(1)-X(2):REM !!!
9020 RETURN