復(fù)合形法源程序
副標(biāo)題：
作者：johhan 文章來源：本站原創(chuàng) 點(diǎn)擊數(shù)：990 更新時(shí)間：2005-8-6

請(qǐng)大家笑納

!**********************************************************************************************************
!* 說明: 1.本程序?yàn)閺?fù)合形法                                                                          *
!*                                                                                                        *
!*       2.程序功能是求解約束最優(yōu)化問題                                                              *
!*   max F(x1,x2,…，xn)                                                               *
!*        s.t.        Gi≤xi≤Hi                                                              *
!*      其中 x1,x2,…，xn為獨(dú)立自變量，                                                              *
!*       xn+1~xm為隱式變量，是x1,x2,…，xn的函數(shù)，Gi，Hi為下界和上界,它們可以是常數(shù)(顯式約束),*
!*  也可以是自變量的函數(shù)(隱式約束)                                                            *
!* 主程序中提供自變量的初始值，輸入已知參數(shù)及打印最后結(jié)果                                            *
!*                                                                                                        *
!*       3.輸入變量說明                                                                              *
!*  N--顯式自變量數(shù)                                                                           *
!*  M--約束組數(shù)                                                                               *
!*  K--構(gòu)成復(fù)合形的頂點(diǎn)數(shù),常用N+1,可多取                                                      *
!*  ITMAX--允許最多迭代次數(shù)                                                                   *
!*  IPRINT--打印控制參數(shù),IPRINT=1,打印中間結(jié)果,IPRINT=0,不打印中間結(jié)果                        *
!*  ALPHA--反射因子,常用1.3                                                                   *
!*  BETA--收斂參數(shù),例如函數(shù)的幅值乘1E-4                                                       *
!*  GAMMA--收斂參數(shù),整數(shù),常用值為5                                                            *
!*  DELTA--顯式約束違反校正,小正數(shù),如X向量幅值乘1E-4                                          *
!*  X(1,J)--自變量初始可行點(diǎn),=1,N                                                             *
!*                                                                                                      *
!*      4.輸出變量說明                                                                               *
!*  F--目標(biāo)函數(shù)最大值勤                                                                       *
!*  X(I)--自變量最優(yōu)值,I=1,N                                                                  *
!*                                                                                                      *
!*      5.使用方法                                                                                   *
!*  1)用戶按照待解問題修改主程序PARAMETER的N,M,K值                                            *
!*  2)在子程序FUNC中給定目標(biāo)函數(shù)                                                              *
!*  3)在子程序CONST中給定Hi和Gi,顯式約束必須放在隱式約束前面                                  *
!*  4)編寫一個(gè)數(shù)據(jù)文件COMPDAT,文件中的數(shù)據(jù)依次為下列變量的值,每個(gè)數(shù)之間用逗號(hào)分開             *  
!*      ITMAX,IPRINT,ALPHA,BETA,GAMMA,DELTA,X(1,1),X(1,2),…，X(1,N)                          *
!*      6.程序內(nèi)容                                                                                   *
!*  1)本程序由一個(gè)主程序和六個(gè)子程序組成，主程序首先給定N,M,K值，然后定義數(shù)組維數(shù)             *
!*   PARAMETER(N=,M=,K=)                                                               *
!*   DIMENSION X(K,M),R(K,N),F(K),G(M),H(M),XC(N)                                      *
!*  主程序中提供自變量的初始值，輸入已知參數(shù)及打印最后結(jié)果面                                  *
!*  2)各子程序的作用分別為了                                                                  *  
!*      CONSX—這是一個(gè)主要子程序，調(diào)用其它于程序及輸出中間結(jié)果                               *
!*      CHECK—檢查所有的點(diǎn)是否滿足約束條件,對(duì)違背約束的點(diǎn)進(jìn)行校正                            *
!*      CENTR—計(jì)算中心點(diǎn)                                                                     *
!*      FUNC —目標(biāo)函數(shù)，由用戶提供                                                           *
!*      CONST—規(guī)定顯式和隱式約束                                                             *
!*      RANDU—產(chǎn)生隨機(jī)數(shù)                                                                     *
!**********************************************************************************************************
 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

!*********************************************************************************************
!* 這是一個(gè)主要子程序，調(diào)用其它于程序及輸出中間結(jié)果                                     *
!*********************************************************************************************
 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

!***********************************************************************************
!* 檢查所有的點(diǎn)是否滿足約束條件,對(duì)違背約束的點(diǎn)進(jìn)行校正。                      *
!***********************************************************************************
 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

!***********************************************************************************
!* 計(jì)算中心點(diǎn)                                                                 *
!***********************************************************************************
 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


!**********************************************************************************
!* 目標(biāo)函數(shù)，由用戶提供                                                      *
!**********************************************************************************
 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


!**********************************************************************************
!* 規(guī)定顯式和隱式約束                                                        *
!**********************************************************************************
 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


!***********************************************************************************
!* 產(chǎn)生隨機(jī)數(shù)                                                                 *
!***********************************************************************************
 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