您所在位置：主题：求助：程序运行到一半出现数学错误
ricky06
[专家分：0] 发布于 2010-06-05 10:50:00
c 复合形法
c N-显式自变量数
c M-约束组数
c K-构成复合形的顶点数，常用N+1，可多取
c ITMAX-允许最大的迭代次数
c IPRINT-打印控制参数，IPRINT=1，打印中间结果
c IPRINT=0,不打印中间结果
c ALPHA-反射因子，常取1.3
c BETA-收敛参数
c GAMMA-收敛参数，整数，常用值5
c DELTA-显式约束违反校正，小正数
c X(1,J)-自变量初始可行点，J=1,N
c F-目标函数最大值
c X(I)-自变量最优值I=1,N
c 程序运行开始，必须回答屏幕提问UNIT2（键入COMPL.DAT),这样，输出的结果将保存在文件COMPL.DAT
C 目标函数在子程序FUNC中由用户 给定
PROGRAM COMPLEXOPTIMIZATIOM
PARAMETER(N=3,M=4,K=6)
DIMENSION X(K,M),R(K,N),F(K),G(M),H(M),XC(N)
INTEGER GAMMA
ITMAX=200
IPRINT=1
ALPHA=1.3
BETA=0.1
GAMMA=5
DELTA=0.0001
DO 110 J=1,N
110 X(1,J)=1.0
IX=2097151
YFL=0.
DO 100 II=2,K
DO 100 JJ=1,N
CALL RANDU(IX,YFL)
R(II,JJ)=YFL
100 CONTINUE
WRITE(2,10)
10 FORMAT(//,18X,'COMPLEX PROGRAM OF BOX')
WRITE(2,11)N,M,K,ITMAX,IPRINT,ALPHA,BETA,GAMMA,DELTA
11 FORMAT(//,2X,'N=',I2,3X,'M=',I2,3X,'K=',I2,3X,'ITMAX=',I4,/,2X,
^ 'IPRINT=',I2,2X,'ALPHA=',F10.4,5X,'BETA=',F10.5,3X,/,2X,'GAMMA=',
^I2,3X,'DELTA=',F10.5)
IF(IPRINT)40,50,40
40 WRITE(2,12)
12 FORMAT(//,2X,'RANDOM NUMBERS')
DO 200 J=2,K
WRITE(2,13)(J,I,R(J,I),I=1,N)
13 FORMAT(/,3(2X,'R(',I2,',',I2,')=',F6.4,2X))
200 CONTINUE
50 CALL CONSX(N,M,K,ITMAX,ALPHA,BETA,GAMMA,DELTA,X,R,F,IT,
^IEV2,G,H,XC,IPRINT)
IF(IT-ITMAX)20,20,30
20 WRITE(2,14)F(IEV2)
14 FORMAT(/,2X,'FINAL VALUE OF THE FUNCTION=',E16.8)
WRITE(2,15)
15 FORMAT(/,2X,'FINAL X VALUES')
DO 300 J=1,N
WRITE(2,16)J,X(IEV2,J)
16 FORMAT(/,2X,'X(',I2,')=',E20.8)
300 CONTINUE
GOTO 999
30 WRITE(2,17)ITMAX
17 FORMAT(/,2X,'THE NUMBER OF ITERATIONS HAS EXCEEDED',I4,10X,
^ 'PROGRAM TERMINATED')
999 STOP
END
!*********************************************************************************************
!* 这是一个主要子程序，调用其它于程序及输出中间结果 *
!*********************************************************************************************
SUBROUTINE CONSX(N,M,K,ITMAX,ALPHA,BETA,GAMMA,DELTA,X,R,F,IT,IEV2,
^G,H,XC,IPRINT)
DIMENSION X(K,M),R(K,N),F(K),G(M),H(M),XC(N)
INTEGER GAMMA
IT=1
KODE=0
IF(M-N)20,20,10
10 KODE=1
20 CONTINUE
DO 40 II=2,K
DO 30 J=1,N
30 X(II,J)=0.0
40 CONTINUE
DO 65 II=2,K
DO 50 J=1,N
I=II
CALL CONST(N,M,K,X,G,H,I)
X(II,J)=G(J)+R(II,J)*(H(J)-G(J))
50 CONTINUE
K1=II
CALL CHECK(N,M,K,X,G,H,I,KODE,XC,DELTA,K1)
IF(II-2)51,51,55
51 IF(IPRINT)52,65,52
52 WRITE(2,18)
18 FORMAT(/,2X,'COORDINATES OF INITIAL COMPLEX')
I0=1
WRITE(2,19)(I0,J,X(I0,J),J=1,N)
19 FORMAT(/,3(2X,'X(',I2,',',I2,')=',1PE13.6))
55 IF(IPRINT)56,65,56
56 WRITE(2,19)(II,J,X(II,J),J=1,N)
65 CONTINUE
K1=K
DO 70 I=1,K
CALL FUNC(N,M,K,X,F,I)
70 CONTINUE
KOUNT=1
IA=0
IF(IPRINT)72,80,72
72 WRITE(2,21)
21 FORMAT(/,2X,'VALUES OF THE FUNCTION')
WRITE(2,22)(J,F(J),J=1,K)
22 FORMAT(/,3(2X,'F(',I2,')=',E13.6))
80 IEV1=1
DO 100 ICM=2,K
IF(F(IEV1)-F(ICM))100,100,90
90 IEV1=ICM
100 CONTINUE
IEV2=1
DO 120 ICM=2,K
IF(F(IEV2)-F(ICM))110,110,120
110 IEV2=ICM
120 CONTINUE
IF(F(IEV2)-(F(IEV1)+BETA))140,130,130
130 KOUNT=1
GOTO 150
140 KOUNT=KOUNT+1
IF(KOUNT-GAMMA)150,240,240
! REPLACEMENT POINT WITH LOWEST FUNCTION VALUE
150 CALL CENTR(N,M,K,IEV1,I,XC,X,K1)
DO 160 JJ=1,N
160 X(IEV1,JJ)=(1.+ALPHA)*(XC(JJ))-ALPHA*(X(IEV1,JJ))
I=IEV1
CALL CHECK(N,M,K,X,G,H,I,KODE,XC,DELTA,K1)
CALL FUNC(N,M,K,X,F,I)
! REPLACEMENT NEW POINT IF IT REPEATS AS LOWEST FUNCTION VALUE
170 IEV2=1
DO 190 ICM=2,K
IF(F(IEV2)-F(ICM))190,190,180
180 IEV2=ICM
190 CONTINUE
IF(IEV2-IEV1)220,200,220
200 DO 210 JJ=1,N
X(IEV1,JJ)=(X(IEV1,JJ)+XC(JJ))/2.
210 CONTINUE
I=IEV1
CALL CHECK(N,M,K,X,G,H,I,KODE,XC,DELTA,K1)
CALL FUNC(N,M,K,X,F,I)
GOTO 170
220 CONTINUE
IF(IPRINT)230,228,230
230 WRITE(2,23)IT
23 FORMAT(//,2X,'ITERATION NUMBER',I5)
WRITE(2,24)
24 FORMAT(/,2X,'COORDINATES OF CORRECTED POINT')
WRITE(2,19)(IEV1,JC,X(IEV1,JC),JC=1,N)
WRITE(2,21)
WRITE(2,22)(I,F(I),I=1,K)
WRITE(2,25)
25 FORMAT(/,2X,'COORDINATES OF CCENTROID')
WRITE(2,26)(JC,XC(JC),JC=1,N)
26 FORMAT(/,3(2X,'X(',I2,',C)=',E14.6,4X))
228 IT=IT+1
IF(IT-ITMAX)80,80,240
240 RETURN
END
!***********************************************************************************
!* 检查所有的点是否满足约束条件,对违背约束的点进行校正。 *
!***********************************************************************************
SUBROUTINE CHECK(N,M,K,X,G,H,I,KODE,XC,DELTA,K1)
! ARGUMENT LIST
! ALL ARGUMENTS DEFINE IN MAIN LINE AND CONSX
DIMENSION X(K,M),G(M),H(M),XC(N)
10 KT=0
CALL CONST(N,M,K,X,G,H,I)
! CHECK AGAINST EXPLICIT CONSTRAINTS
DO 50 J=1,N
IF(X(I,J)-G(J))20,20,30
20 X(I,J)=G(J)+DELTA
GOTO 50
30 IF(H(J)-X(I,J))40,40,50
40 X(I,J)=H(J)-DELTA
50 CONTINUE
IF(KODE)110,110,60
! CHECK AGAINST THE IMPLICIT CONSTTRAINTS
60 NN=N+1
DO 100 J=NN,M
CALL CONST(N,M,K,X,G,H,I)
IF(X(I,J)-G(J))80,70,70
70 IF(H(J)-X(I,J))80,100,100
80 IEV1=I
KT=1
CALL CENTR(N,M,K,IEV1,I,XC,X,K1)
DO 90 JJ=1,N
X(I,JJ)=(X(I,JJ)+XC(JJ))/2
90 CONTINUE
100 CONTINUE
IF(KT)110,110,10
110 RETURN
END
!***********************************************************************************
!* 计算中心点 *
!***********************************************************************************
SUBROUTINE CENTR(N,M,K,IEV1,I,XC,X,K1)
DIMENSION X(K,M),XC(N)
DO 20 J=1,N
XC(J)=0.
DO 10 IL=1,K1
10 XC(J)=XC(J)+X(IL,J)
RK=K1
20 XC(J)=(XC(J)-X(IEV1,J))/(RK-1.)
RETURN
END
!**********************************************************************************
!* 目标函数，由用户提供 *
!**********************************************************************************
SUBROUTINE FUNC(N,M,K,X,F,I)
DIMENSION X(K,M),F(K)
! OBJECTIVE FUNCTION
F(I)=(2*8*(
^0.5*(X(I,1)+X(I,3))*SQRT((94-2*X(I,3))*X(I,1)+47*47-X(I,3)*X(I,3))
^+0.5*(X(I,2)+X(I,3))*SQRT((54-2*X(I,3))*X(I,2)+27*27-X(I,3)*X(I,3)
^)-0.5*27*27*ASIN((X(I,2)+X(I,3))/(X(I,2)+27))
^-0.5*47*47*ASIN((X(I,1)+X(I,3))/(X(I,1)+47))
^-0.5*X(I,2)*X(I,2)*ACOS((X(I,2)+X(I,3))/(X(I,2)+27))
^-0.5*X(I,1)*X(I,1)*ACOS((X(I,1)+X(I,3))/(X(I,1)+47))
^+X(I,3)*(105-SQRT((94-2*X(I,3))*X(I,1)+47*47-X(I,3)*X(I,3))
^-SQRT((54-2*X(I,3))*X(I,1)+27*27-X(I,3)*X(I,3))))
^-10*(X(I,3)-3))
RETURN
END
!**********************************************************************************
!* 规定显式和隐式约束 *
!**********************************************************************************
SUBROUTINE CONST(N,M,K,X,G,H,I)
DIMENSION X(K,M),G(M),H(M)
! CONSSTRAAINTS LIMITS ND FUNCTION
G(1)=39.
G(2)=19.
G(3)=3.
G(4)=0.
H(1)=100.
H(2)=100.
H(3)=10.
H(4)=523.
X(I,4)=12517.2/(2*(24+3*((X(I,3-3)))
RETURN
END
!***********************************************************************************
!* 产生随机数 *
!***********************************************************************************
SUBROUTINE RANDU(IX,YFL)
IF(YFL.NE.0.0)GOTO 1
IM=2**21
IC=2**10-3
AX=FLOAT(IX)
AM=FLOAT(IM)
AC=FLOAT(IC)
YFL=AX/AM
1 YFL=AC*YFL
YFL=YFL-FLOAT(IFIX(YFL))
RETURN
END
c N-显式自变量数
c M-约束组数
c K-构成复合形的顶点数，常用N+1，可多取
c ITMAX-允许最大的迭代次数
c IPRINT-打印控制参数，IPRINT=1，打印中间结果
c IPRINT=0,不打印中间结果
c ALPHA-反射因子，常取1.3
c BETA-收敛参数
c GAMMA-收敛参数，整数，常用值5
c DELTA-显式约束违反校正，小正数
c X(1,J)-自变量初始可行点，J=1,N
c F-目标函数最大值
c X(I)-自变量最优值I=1,N
c 程序运行开始，必须回答屏幕提问UNIT2（键入COMPL.DAT),这样，输出的结果将保存在文件COMPL.DAT
C 目标函数在子程序FUNC中由用户 给定
PROGRAM COMPLEXOPTIMIZATIOM
PARAMETER(N=3,M=4,K=6)
DIMENSION X(K,M),R(K,N),F(K),G(M),H(M),XC(N)
INTEGER GAMMA
ITMAX=200
IPRINT=1
ALPHA=1.3
BETA=0.1
GAMMA=5
DELTA=0.0001
DO 110 J=1,N
110 X(1,J)=1.0
IX=2097151
YFL=0.
DO 100 II=2,K
DO 100 JJ=1,N
CALL RANDU(IX,YFL)
R(II,JJ)=YFL
100 CONTINUE
WRITE(2,10)
10 FORMAT(//,18X,'COMPLEX PROGRAM OF BOX')
WRITE(2,11)N,M,K,ITMAX,IPRINT,ALPHA,BETA,GAMMA,DELTA
11 FORMAT(//,2X,'N=',I2,3X,'M=',I2,3X,'K=',I2,3X,'ITMAX=',I4,/,2X,
^ 'IPRINT=',I2,2X,'ALPHA=',F10.4,5X,'BETA=',F10.5,3X,/,2X,'GAMMA=',
^I2,3X,'DELTA=',F10.5)
IF(IPRINT)40,50,40
40 WRITE(2,12)
12 FORMAT(//,2X,'RANDOM NUMBERS')
DO 200 J=2,K
WRITE(2,13)(J,I,R(J,I),I=1,N)
13 FORMAT(/,3(2X,'R(',I2,',',I2,')=',F6.4,2X))
200 CONTINUE
50 CALL CONSX(N,M,K,ITMAX,ALPHA,BETA,GAMMA,DELTA,X,R,F,IT,
^IEV2,G,H,XC,IPRINT)
IF(IT-ITMAX)20,20,30
20 WRITE(2,14)F(IEV2)
14 FORMAT(/,2X,'FINAL VALUE OF THE FUNCTION=',E16.8)
WRITE(2,15)
15 FORMAT(/,2X,'FINAL X VALUES')
DO 300 J=1,N
WRITE(2,16)J,X(IEV2,J)
16 FORMAT(/,2X,'X(',I2,')=',E20.8)
300 CONTINUE
GOTO 999
30 WRITE(2,17)ITMAX
17 FORMAT(/,2X,'THE NUMBER OF ITERATIONS HAS EXCEEDED',I4,10X,
^ 'PROGRAM TERMINATED')
999 STOP
END
!*********************************************************************************************
!* 这是一个主要子程序，调用其它于程序及输出中间结果 *
!*********************************************************************************************
SUBROUTINE CONSX(N,M,K,ITMAX,ALPHA,BETA,GAMMA,DELTA,X,R,F,IT,IEV2,
^G,H,XC,IPRINT)
DIMENSION X(K,M),R(K,N),F(K),G(M),H(M),XC(N)
INTEGER GAMMA
IT=1
KODE=0
IF(M-N)20,20,10
10 KODE=1
20 CONTINUE
DO 40 II=2,K
DO 30 J=1,N
30 X(II,J)=0.0
40 CONTINUE
DO 65 II=2,K
DO 50 J=1,N
I=II
CALL CONST(N,M,K,X,G,H,I)
X(II,J)=G(J)+R(II,J)*(H(J)-G(J))
50 CONTINUE
K1=II
CALL CHECK(N,M,K,X,G,H,I,KODE,XC,DELTA,K1)
IF(II-2)51,51,55
51 IF(IPRINT)52,65,52
52 WRITE(2,18)
18 FORMAT(/,2X,'COORDINATES OF INITIAL COMPLEX')
I0=1
WRITE(2,19)(I0,J,X(I0,J),J=1,N)
19 FORMAT(/,3(2X,'X(',I2,',',I2,')=',1PE13.6))
55 IF(IPRINT)56,65,56
56 WRITE(2,19)(II,J,X(II,J),J=1,N)
65 CONTINUE
K1=K
DO 70 I=1,K
CALL FUNC(N,M,K,X,F,I)
70 CONTINUE
KOUNT=1
IA=0
IF(IPRINT)72,80,72
72 WRITE(2,21)
21 FORMAT(/,2X,'VALUES OF THE FUNCTION')
WRITE(2,22)(J,F(J),J=1,K)
22 FORMAT(/,3(2X,'F(',I2,')=',E13.6))
80 IEV1=1
DO 100 ICM=2,K
IF(F(IEV1)-F(ICM))100,100,90
90 IEV1=ICM
100 CONTINUE
IEV2=1
DO 120 ICM=2,K
IF(F(IEV2)-F(ICM))110,110,120
110 IEV2=ICM
120 CONTINUE
IF(F(IEV2)-(F(IEV1)+BETA))140,130,130
130 KOUNT=1
GOTO 150
140 KOUNT=KOUNT+1
IF(KOUNT-GAMMA)150,240,240
! REPLACEMENT POINT WITH LOWEST FUNCTION VALUE
150 CALL CENTR(N,M,K,IEV1,I,XC,X,K1)
DO 160 JJ=1,N
160 X(IEV1,JJ)=(1.+ALPHA)*(XC(JJ))-ALPHA*(X(IEV1,JJ))
I=IEV1
CALL CHECK(N,M,K,X,G,H,I,KODE,XC,DELTA,K1)
CALL FUNC(N,M,K,X,F,I)
! REPLACEMENT NEW POINT IF IT REPEATS AS LOWEST FUNCTION VALUE
170 IEV2=1
DO 190 ICM=2,K
IF(F(IEV2)-F(ICM))190,190,180
180 IEV2=ICM
190 CONTINUE
IF(IEV2-IEV1)220,200,220
200 DO 210 JJ=1,N
X(IEV1,JJ)=(X(IEV1,JJ)+XC(JJ))/2.
210 CONTINUE
I=IEV1
CALL CHECK(N,M,K,X,G,H,I,KODE,XC,DELTA,K1)
CALL FUNC(N,M,K,X,F,I)
GOTO 170
220 CONTINUE
IF(IPRINT)230,228,230
230 WRITE(2,23)IT
23 FORMAT(//,2X,'ITERATION NUMBER',I5)
WRITE(2,24)
24 FORMAT(/,2X,'COORDINATES OF CORRECTED POINT')
WRITE(2,19)(IEV1,JC,X(IEV1,JC),JC=1,N)
WRITE(2,21)
WRITE(2,22)(I,F(I),I=1,K)
WRITE(2,25)
25 FORMAT(/,2X,'COORDINATES OF CCENTROID')
WRITE(2,26)(JC,XC(JC),JC=1,N)
26 FORMAT(/,3(2X,'X(',I2,',C)=',E14.6,4X))
228 IT=IT+1
IF(IT-ITMAX)80,80,240
240 RETURN
END
!***********************************************************************************
!* 检查所有的点是否满足约束条件,对违背约束的点进行校正。 *
!***********************************************************************************
SUBROUTINE CHECK(N,M,K,X,G,H,I,KODE,XC,DELTA,K1)
! ARGUMENT LIST
! ALL ARGUMENTS DEFINE IN MAIN LINE AND CONSX
DIMENSION X(K,M),G(M),H(M),XC(N)
10 KT=0
CALL CONST(N,M,K,X,G,H,I)
! CHECK AGAINST EXPLICIT CONSTRAINTS
DO 50 J=1,N
IF(X(I,J)-G(J))20,20,30
20 X(I,J)=G(J)+DELTA
GOTO 50
30 IF(H(J)-X(I,J))40,40,50
40 X(I,J)=H(J)-DELTA
50 CONTINUE
IF(KODE)110,110,60
! CHECK AGAINST THE IMPLICIT CONSTTRAINTS
60 NN=N+1
DO 100 J=NN,M
CALL CONST(N,M,K,X,G,H,I)
IF(X(I,J)-G(J))80,70,70
70 IF(H(J)-X(I,J))80,100,100
80 IEV1=I
KT=1
CALL CENTR(N,M,K,IEV1,I,XC,X,K1)
DO 90 JJ=1,N
X(I,JJ)=(X(I,JJ)+XC(JJ))/2
90 CONTINUE
100 CONTINUE
IF(KT)110,110,10
110 RETURN
END
!***********************************************************************************
!* 计算中心点 *
!***********************************************************************************
SUBROUTINE CENTR(N,M,K,IEV1,I,XC,X,K1)
DIMENSION X(K,M),XC(N)
DO 20 J=1,N
XC(J)=0.
DO 10 IL=1,K1
10 XC(J)=XC(J)+X(IL,J)
RK=K1
20 XC(J)=(XC(J)-X(IEV1,J))/(RK-1.)
RETURN
END
!**********************************************************************************
!* 目标函数，由用户提供 *
!**********************************************************************************
SUBROUTINE FUNC(N,M,K,X,F,I)
DIMENSION X(K,M),F(K)
! OBJECTIVE FUNCTION
F(I)=(2*8*(
^0.5*(X(I,1)+X(I,3))*SQRT((94-2*X(I,3))*X(I,1)+47*47-X(I,3)*X(I,3))
^+0.5*(X(I,2)+X(I,3))*SQRT((54-2*X(I,3))*X(I,2)+27*27-X(I,3)*X(I,3)
^)-0.5*27*27*ASIN((X(I,2)+X(I,3))/(X(I,2)+27))
^-0.5*47*47*ASIN((X(I,1)+X(I,3))/(X(I,1)+47))
^-0.5*X(I,2)*X(I,2)*ACOS((X(I,2)+X(I,3))/(X(I,2)+27))
^-0.5*X(I,1)*X(I,1)*ACOS((X(I,1)+X(I,3))/(X(I,1)+47))
^+X(I,3)*(105-SQRT((94-2*X(I,3))*X(I,1)+47*47-X(I,3)*X(I,3))
^-SQRT((54-2*X(I,3))*X(I,1)+27*27-X(I,3)*X(I,3))))
^-10*(X(I,3)-3))
RETURN
END
!**********************************************************************************
!* 规定显式和隐式约束 *
!**********************************************************************************
SUBROUTINE CONST(N,M,K,X,G,H,I)
DIMENSION X(K,M),G(M),H(M)
! CONSSTRAAINTS LIMITS ND FUNCTION
G(1)=39.
G(2)=19.
G(3)=3.
G(4)=0.
H(1)=100.
H(2)=100.
H(3)=10.
H(4)=523.
X(I,4)=12517.2/(2*(24+3*((X(I,3-3)))
RETURN
END
!***********************************************************************************
!* 产生随机数 *
!***********************************************************************************
SUBROUTINE RANDU(IX,YFL)
IF(YFL.NE.0.0)GOTO 1
IM=2**21
IC=2**10-3
AX=FLOAT(IX)
AM=FLOAT(IM)
AC=FLOAT(IC)
YFL=AX/AM
1 YFL=AC*YFL
YFL=YFL-FLOAT(IFIX(YFL))
RETURN
END