fuhexingfa.txt

复合形法源程序

!* 说明: 1.本程序为复合形法                                                                          *
!*                                                                                                        *
!*       2.程序功能是求解约束最优化问题                                                              *
!*   max F(x1,x2,…，xn)                                                               *
!*        s.t.        Gi≤xi≤Hi                                                              *
!*      其中 x1,x2,…，xn为独立自变量，                                                              *
!*       xn+1~xm为隐式变量，是x1,x2,…，xn的函数，Gi，Hi为下界和上界,它们可以是常数(显式约束),*
!*  也可以是自变量的函数(隐式约束)                                                            *
!* 主程序中提供自变量的初始值，输入已知参数及打印最后结果                                            *
!*                                                                                                        *
!*       3.输入变量说明                                                                              *
!*  N--显式自变量数                                                                           *
!*  M--约束组数                                                                               *
!*  K--构成复合形的顶点数,常用N+1,可多取                                                      *
!*  ITMAX--允许最多迭代次数                                                                   *
!*  IPRINT--打印控制参数,IPRINT=1,打印中间结果,IPRINT=0,不打印中间结果                        *
!*  ALPHA--反射因子,常用1.3                                                                   *
!*  BETA--收敛参数,例如函数的幅值乘1E-4                                                       *
!*  GAMMA--收敛参数,整数,常用值为5                                                            *
!*  DELTA--显式约束违反校正,小正数,如X向量幅值乘1E-4                                          *
!*  X(1,J)--自变量初始可行点,=1,N                                                             *
!*                                                                                                      *
!*      4.输出变量说明                                                                               *
!*  F--目标函数最大值勤                                                                       *
!*  X(I)--自变量最优值,I=1,N                                                                  *
!*                                                                                                      *
!*      5.使用方法                                                                                   *
!*  1)用户按照待解问题修改主程序PARAMETER的N,M,K值                                            *
!*  2)在子程序FUNC中给定目标函数                                                              *
!*  3)在子程序CONST中给定Hi和Gi,显式约束必须放在隐式约束前面                                  *
!*  4)编写一个数据文件COMPDAT,文件中的数据依次为下列变量的值,每个数之间用逗号分开             *   
!*      ITMAX,IPRINT,ALPHA,BETA,GAMMA,DELTA,X(1,1),X(1,2),…，X(1,N)                          *
!*      6.程序内容                                                                                   *
!*  1)本程序由一个主程序和六个子程序组成，主程序首先给定N,M,K值，然后定义数组维数             *
!*   PARAMETER(N=,M=,K=)                                                               *
!*   DIMENSION X(K,M),R(K,N),F(K),G(M),H(M),XC(N)                                      *
!*  主程序中提供自变量的初始值，输入已知参数及打印最后结果面                                  *
!*  2)各子程序的作用分别为了                                                                  *   
!*      CONSX—这是一个主要子程序，调用其它于程序及输出中间结果                               *
!*      CHECK—检查所有的点是否满足约束条件,对违背约束的点进行校正                            *
!*      CENTR—计算中心点                                                                     *
!*      FUNC —目标函数，由用户提供                                                           *
!*      CONST—规定显式和隐式约束                                                             *
!*      RANDU—产生随机数                                                                     *
!**********************************************************************************************************
 PROGRAM COMPLEX
 PARAMETER(N=3,M=4,K=6)
 DIMENSION X(K,M),R(K,N),F(K),G(M),H(M),XC(N)
 INTEGER GAMMA
 OPEN(4,FILE='COMPDATA')
 READ(4,*)ITMAX,IPRINT,ALPHA,BETA,GAMMA,DELTA
 READ(4,*)(X(1,J),J=1,N)
 IX=2097151
 YFL=0.
 DO 100 II=2,K
 DO 100 JJ=1,N
 CALL RANDU(IX,YFL)
 R(II,JJ)=YFL
100 CONTINUE
 OPEN(2,FILE='COMPDAT.OUT',STATUS='UNKNOWN')
 WRITE(2,10)
10 FORMAT(//,10X,'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,&
  /,2X,'GAMMA=',I2,2X,'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)=X(I,1)*X(I,2)*X(I,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)=0.
 G(2)=0.
 G(3)=0.
 G(4)=0.
 H(1)=42.
 H(2)=42.
 H(3)=42.
 H(4)=72.
 X(I,4)=X(I,1)+2.*X(I,2)+2*X(I,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