基本原理
SCE-UA算法是Qingyun Duan（段青云）、Soroosh Sorooshian 和Vijai Gupta等开发的一个具有复合优化策略的优化算法（Duan等，1992）。具体原理可以参考文献[1-3]，以下图片摘录自王书功博士在《水文模型参数估计方法及参数估计不确定性研究》2.2.2.3章节（P101-108）中的介绍[5]。这里的CCE算法步骤分解不清楚，应当参考文献[2]。
参考文献
[1]Duan, Q., Sorooshian, S., and Gupta, V. (1992), Effective and efficient global optimization for conceptual rainfall-runoff models, Water Resour. Res., 28( 4), 1015– 1031, doi:10.1029/91WR02985.
[2]Duan, Q.Y., Gupta, V.K. & Sorooshian, S. Shuffled complex evolution approach for effective and efficient global minimization. J Optim Theory Appl 76, 501–521 (1993). https://doi.org/10.1007/BF00939380
[3]Duan, Q., Sorooshian, S., & Gupta, V. K. (1994). Optimal use of the SCE-UA global optimization method for calibrating watershed models. Journal of Hydrology, 158(3-4), 265-284. https://doi.org/10.1016/0022-1694(94)90057-4
[4]王书功. 水文模型参数估计方法及参数估计不确定性研究[M]. 河南:黄河水利出版社,2010.
[5]王书功. 水文模型参数估计方法及参数估计不确定性研究[D]. 北京:中国科学院研究生院,2006.
原始代码
以下代码仅供参考，都无法成功运行。
可以成功运行的C++代码，请从我的另一篇博客《【算法】03 SCE-UA算法C++实现》下载
博客作者简介
很高兴认识您！
我叫卢家波，河海大学水文学及水资源博士研究生，研究兴趣为高效洪水淹没预测、洪水灾害预警、机器学习、替代模型和降阶模型。
变化环境下，极端洪水事件多发，我希望能通过研究为水灾害防御做出贡献，为人民服务。
欢迎交流讨论和研究合作，vx Jiabo_Lu。
主页 https://lujiabo98.github.io
简历 https://lujiabo98.github.io/file/CV_JiaboLu_zh.pdf
博客 https://blog.csdn.net/weixin_43012724?type=blog
来信请说明博客标题及链接，谢谢。
源码
sceua.for：
SUBROUTINE SCEUA(A,AF,BL,BU,NOPT,MAXN,KSTOP,PCENTO,ISEED,
& NGS,NPG,NPS,NSPL,MINGS,INIFLG,IPRINT)
C
C
C SHUFFLED COMPLEX EVOLUTION METHOD FOR GLOBAL OPTIMIZATION
C -- VERSION 2.1
C
C BY QINGYUN DUAN
C DEPARTMENT OF HYDROLOGY & WATER RESOURCES
C UNIVERSITY OF ARIZONA, TUCSON, AZ 85721
C (602) 621-9360, EMAIL: DUAN@HWR.ARIZONA.EDU
C
C WRITTEN IN OCTOBER 1990.
C REVISED IN AUGUST 1991
C REVISED IN APRIL 1992
C
C STATEMENT BY AUTHOR:
C --------------------
C
C THIS GENERAL PURPOSE GLOBAL OPTIMIZATION PROGRAM IS DEVELOPED AT
C THE DEPARTMENT OF HYDROLOGY & WATER RESOURCES OF THE UNIVERSITY
C OF ARIZONA. FURTHER INFORMATION REGARDING THE SCE-UA METHOD CAN
C BE OBTAINED FROM DR. Q. DUAN, DR. S. SOROOSHIAN OR DR. V.K. GUPTA
C AT THE ADDRESS AND PHONE NUMBER LISTED ABOVE. WE REQUEST ALL
C USERS OF THIS PROGRAM MAKE PROPER REFERENCE TO THE PAPER ENTITLED
C 'EFFECTIVE AND EFFICIENT GLOBAL OPTIMIZATION FOR CONCEPTUAL
C RAINFALL-RUNOFF MODELS' BY DUAN, Q., S. SOROOSHIAN, AND V.K. GUPTA,
C WATER RESOURCES RESEARCH, VOL 28(4), PP.1015-1031, 1992.
C
C
C LIST OF INPUT ARGUEMENT VARIABLES
C
C A(.) = INITIAL PARAMETER SET
C BL(.) = LOWER BOUND ON PARAMETERS
C BU(.) = UPPER BOUND ON PARAMETERS
C NOPT = NUMBER OF PARAMETERS TO BE OPTIMIZED
C
C
C LIST OF SCE ALGORITHMIC CONTROL PARAMETERS:
C
C NGS = NUMBER OF COMPLEXES IN THE INITIAL POPULATION
C NPG = NUMBER OF POINTS IN EACH COMPLEX
C NPT = TOTAL NUMBER OF POINTS IN INITIAL POPULATION (NPT=NGS*NPG)
C NPS = NUMBER OF POINTS IN A SUB-COMPLEX
C NSPL = NUMBER OF EVOLUTION STEPS ALLOWED FOR EACH COMPLEX BEFORE
C COMPLEX SHUFFLING
C MINGS = MINIMUM NUMBER OF COMPLEXES REQUIRED, IF THE NUMBER OF
C COMPLEXES IS ALLOWED TO REDUCE AS THE OPTIMIZATION PROCEEDS
C ISEED = INITIAL RANDOM SEED
C INIFLG = FLAG ON WHETHER TO INCLUDE THE INITIAL POINT IN POPULATION
C = 0, NOT INCLUDED
C = 1, INCLUDED
C IPRINT = FLAG FOR CONTROLLING PRINT-OUT AFTER EACH SHUFFLING LOOP
C = 0, PRINT INFORMATION ON THE BEST POINT OF THE POPULATION
C = 1, PRINT INFORMATION ON EVERY POINT OF THE POPULATION
C
C
C CONVERGENCE CHECK PARAMETERS
C
C MAXN = MAX NO. OF TRIALS ALLOWED BEFORE OPTIMIZATION IS TERMINATED
C KSTOP = NUMBER OF SHUFFLING LOOPS IN WHICH THE CRITERION VALUE MUST
C CHANG BY THE GIVEN PERCENTAGE BEFORE OPTIMIZATION IS TERMINATED
C PCENTO = PERCENTAGE BY WHICH THE CRITERION VALUE MUST CHANGE IN
C GIVEN NUMBER OF SHUFFLING LOOPS
C IPCNVG = FLAG INDICATING WHETHER PARAMETER CONVERGENCE IS REACHED
C (I.E., CHECK IF GNRNG IS LESS THAN 0.001)
C = 0, PARAMETER CONVERGENCE NOT SATISFIED
C = 1, PARAMETER CONVERGENCE SATISFIED
C
C
C LIST OF LOCAL VARIABLES
C X(.,.) = COORDINATES OF POINTS IN THE POPULATION
C XF(.) = FUNCTION VALUES OF X(.,.)
C XX(.) = COORDINATES OF A SINGLE POINT IN X
C CX(.,.) = COORDINATES OF POINTS IN A COMPLEX
C CF(.) = FUNCTION VALUES OF CX(.,.)
C S(.,.) = COORDINATES OF POINTS IN THE CURRENT SIMPLEX
C SF(.) = FUNCTION VALUES OF S(.,.)
C BESTX(.) = BEST POINT AT CURRENT SHUFFLING LOOP
C BESTF = FUNCTION VALUE OF BESTX(.)
C WORSTX(.) = WORST POINT AT CURRENT SHUFFLING LOOP
C WORSTF = FUNCTION VALUE OF WORSTX(.)
C XNSTD(.) = STANDARD DEVIATION OF PARAMETERS IN THE POPULATION
C GNRNG = NORMALIZED GEOMETRIC MEAN OF PARAMETER RANGES
C LCS(.) = INDICES LOCATING POSITION OF S(.,.) IN X(.,.)
C BOUND(.) = BOUND ON ITH VARIABLE BEING OPTIMIZED
C NGS1 = NUMBER OF COMPLEXES IN CURRENT POPULATION
C NGS2 = NUMBER OF COMPLEXES IN LAST POPULATION
C ISEED1 = CURRENT RANDOM SEED
C CRITER(.) = VECTOR CONTAINING THE BEST CRITERION VALUES OF THE LAST
C 10 SHUFFLING LOOPS
C
CHARACTER*4 XNAME(16)
C
C ARRAYS FROM THE INPUT DATA
DIMENSION A(16),BL(16),BU(16)
C
C LOCAL ARRAYS
DIMENSION X(2000,16),XX(16),BESTX(16),WORSTX(16),XF(2000)
DIMENSION S(50,16),SF(50),LCS(50),CX(2000,16),CF(2000)
DIMENSION XNSTD(16),BOUND(16),CRITER(20),UNIT(16)
C
COMMON/IOPARS/ICNTRL,IOUT,IDAT,IWBAL,ISCE,IPE,IPC,IDET
C
DATA XNAME /' X1',' X2',' X3',' X4',' X5',' X6',' X7',
&' X8',' X9',' X10',' X11',' X12',' X13',' X14',' X15',' X16'/
C
C INITIALIZE VARIABLES
WRITE(ISCE,400)
400 FORMAT(//,2X,50(1H=),/,2X,'ENTER THE SHUFFLED COMPLEX EVOLUTION',
& ' GLOBAL SEARCH',/,2X,50(1H=))
WRITE(*,400)
NLOOP = 0
LOOP = 0
IGS = 0
C
C INITIALIZE RANDOM SEED TO A NEGATIVE INTEGER
ISEED1 = -ABS(ISEED)
C
C COMPUTE THE TOTAL NUMBER OF POINTS IN INITIAL POPUALTION
NPT = NGS * NPG
NGS1 = NGS
NPT1 = NPT
C
C COMPUTE THE BOUND FOR PARAMETERS BEING OPTIMIZED
DO J = 1, NOPT
BOUND(J) = BU(J) - BL(J)
UNIT(J) = 1.0
END DO
C
C COMPUTE THE FUNCTION VALUE OF THE INITIAL POINT
FA = FUNCTN(NOPT,A)
C
C PRINT THE INITIAL POINT AND ITS CRITERION VALUE
WRITE(ISCE,500)
WRITE(*, 500)
WRITE(ISCE,510) (XNAME(J),J=1,NOPT)
WRITE(*, 510) (XNAME(J),J=1,NOPT)
WRITE(ISCE,520) FA,(A(J),J=1,NOPT)
WRITE(*, 520) FA,(A(J),J=1,NOPT)
IF (MAXN .EQ. 1) GO TO 10000
C
C GENERATE AN INITIAL SET OF NPT1 POINTS IN THE PARAMETER SPACE
C IF INIFLG IS EQUAL TO 1, SET X(1,.) TO INITIAL POINT A(.)
IF (INIFLG .EQ. 1) THEN
DO J = 1, NOPT
X(1,J) = A(J)
END DO
XF(1) = FA
C
C ELSE, GENERATE A POINT RANDOMLY AND SET IT EQUAL TO X(1,.)
ELSE
CALL GETPNT(NOPT,1,ISEED1,XX,BL,BU,UNIT,BL)
DO J=1, NOPT
X(1,J) = XX(J)
END DO
XF(1) = FUNCTN(NOPT,XX)
END IF
ICALL = 1
IF (ICALL .GE. MAXN) GO TO 9000
C
C GENERATE NPT1-1 RANDOM POINTS DISTRIBUTED UNIFORMLY IN THE PARAMETER
C SPACE, AND COMPUTE THE CORRESPONDING FUNCTION VALUES
DO I = 2, NPT1
CALL GETPNT(NOPT,1,ISEED1,XX,BL,BU,UNIT,BL)
DO J = 1, NOPT
X(I,J) = XX(J)
END DO
XF(I) = FUNCTN(NOPT,XX)
ICALL = ICALL + 1
IF (ICALL .GE. MAXN) THEN
NPT1 = I
GO TO 45
END IF
END DO
C
C ARRANGE THE POINTS IN ORDER OF INCREASING FUNCTION VALUE
45 CALL SORT(NPT1,NOPT,X,XF)
C
C RECORD THE BEST AND WORST POINTS
DO J = 1, NOPT
BESTX(J) = X(1,J)
WORSTX(J) = X(NPT1,J)
END DO
BESTF = XF(1)
WORSTF = XF(NPT1)
C
C COMPUTE THE PARAMETER RANGE FOR THE INITIAL POPULATION
CALL PARSTT(NPT1,NOPT,X,XNSTD,BOUND,GNRNG,IPCNVG)
C
C PRINT THE RESULTS FOR THE INITIAL POPULATION
WRITE(ISCE,600)
WRITE(*, 600)
WRITE(ISCE,610) (XNAME(J),J=1,NOPT)
WRITE(*, 610) (XNAME(J),J=1,NOPT)
WRITE(ISCE,630) NLOOP,ICALL,NGS1,BESTF,WORSTF,GNRNG,
& (BESTX(J),J=1,NOPT)
WRITE(*, 630) NLOOP,ICALL,NGS1,BESTF,WORSTF,GNRNG,
& (BESTX(J),J=1,NOPT)
IF (IPRINT .EQ. 1) THEN
WRITE(ISCE,650) NLOOP
DO I = 1, NPT1
WRITE(ISCE,660) XF(I),GNRNG,(X(I,J),J=1,NOPT)
END DO
END IF
C
IF (ICALL .GE. MAXN) GO TO 9000
IF (IPCNVG .EQ. 1) GO TO 9200
C
C BEGIN THE MAIN LOOP ----------------
1000 CONTINUE
NLOOP = NLOOP + 1
C
C BEGIN LOOP ON COMPLEXES
DO 3000 IGS = 1, NGS1
C
C ASSIGN POINTS INTO COMPLEXES
DO K1 = 1, NPG
K2 = (K1-1) * NGS1 + IGS
DO J = 1, NOPT
CX(K1,J) = X(K2,J)
END DO
CF(K1) = XF(K2)
END DO
C
C BEGIN INNER LOOP - RANDOM SELECTION OF SUB-COMPLEXES ---------------
DO 2000 LOOP = 1, NSPL
C
C CHOOSE A SUB-COMPLEX (NPS POINTS) ACCORDING TO A LINEAR
C PROBABILITY DISTRIBUTION
IF (NPS .EQ. NPG) THEN
DO K = 1, NPS
LCS(K) = K
END DO
GO TO 85
END IF
C
RAND = RAN1(ISEED1)
LCS(1) = 1 + INT(NPG + 0.5 - SQRT( (NPG+.5)**2 -
& NPG * (NPG+1) * RAND ))
DO K = 2, NPS
60 RAND = RAN1(ISEED1)
LPOS = 1 + INT(NPG + 0.5 - SQRT((NPG+.5)**2 -
& NPG * (NPG+1) * RAND ))
DO K1 = 1, K-1
IF (LPOS .EQ. LCS(K1)) GO TO 60
END DO
LCS(K) = LPOS
END DO
C
C ARRANGE THE SUB-COMPLEX IN ORDER OF INCEASING FUNCTION VALUE
CALL SORT1(NPS,LCS)
C
C CREATE THE SUB-COMPLEX ARRAYS
85 DO K = 1, NPS
DO J = 1, NOPT
S(K,J) = CX(LCS(K),J)
END DO
SF(K) = CF(LCS(K))
END DO
C
C USE THE SUB-COMPLEX TO GENERATE NEW POINT(S)
CALL CCE(NOPT,NPS,S,SF,BL,BU,XNSTD,ICALL,MAXN,ISEED1)
C
C IF THE SUB-COMPLEX IS ACCEPTED, REPLACE THE NEW SUB-COMPLEX
C INTO THE COMPLEX
DO K = 1, NPS
DO J = 1, NOPT
CX(LCS(K),J) = S(K,J)
END DO
CF(LCS(K)) = SF(K)
END DO
C
C SORT THE POINTS
CALL SORT(NPG,NOPT,CX,CF)
C
C IF MAXIMUM NUMBER OF RUNS EXCEEDED, BREAK OUT OF THE LOOP
IF (ICALL .GE. MAXN) GO TO 2222
C
C END OF INNER LOOP ------------
2000 CONTINUE
2222 CONTINUE
C
C REPLACE THE NEW COMPLEX INTO ORIGINAL ARRAY X(.,.)
DO K1 = 1, NPG
K2 = (K1-1) * NGS1 + IGS
DO J = 1, NOPT
X(K2,J) = CX(K1,J)
END DO
XF(K2) = CF(K1)
END DO
IF (ICALL .GE. MAXN) GO TO 3333
C
C END LOOP ON COMPLEXES
3000 CONTINUE
C
C RE-SORT THE POINTS
3333 CALL SORT(NPT1,NOPT,X,XF)
C
C RECORD THE BEST AND WORST POINTS
DO J = 1, NOPT
BESTX(J) = X(1,J)
WORSTX(J) = X(NPT1,J)
END DO
BESTF = XF(1)
WORSTF = XF(NPT1)
C
C TEST THE POPULATION FOR PARAMETER CONVERGENCE
CALL PARSTT(NPT1,NOPT,X,XNSTD,BOUND,GNRNG,IPCNVG)
C
C PRINT THE RESULTS FOR CURRENT POPULATION
WRITE(ISCE,630) NLOOP,ICALL,NGS1,BESTF,WORSTF,GNRNG,
& (BESTX(J),J=1,NOPT)
WRITE(*,630) NLOOP,ICALL,NGS1,BESTF,WORSTF,GNRNG,
& (BESTX(J),J=1,NOPT)
IF (IPRINT .EQ. 1) THEN
WRITE(ISCE,650) NLOOP
DO I = 1, NPT1
WRITE(ISCE,660) XF(I),GNRNG,(X(I,J),J=1,NOPT)
END DO
END IF
C
C TEST IF MAXIMUM NUMBER OF FUNCTION EVALUATIONS EXCEEDED
IF (ICALL .GE. MAXN) GO TO 9000
C
C COMPUTE THE COUNT ON SUCCESSIVE LOOPS W/O FUNCTION IMPROVEMENT
CRITER(20) = BESTF
IF (NLOOP .LT. (KSTOP+1)) GO TO 132
DENOMI = ABS(CRITER(20-KSTOP) + CRITER(20)) / 2.
TIMEOU = ABS(CRITER(20-KSTOP) - CRITER(20)) / DENOMI
IF (TIMEOU .LT. PCENTO) GO TO 9100
132 CONTINUE
DO L = 1, 19
CRITER(L) = CRITER(L+1)
END DO
C
C IF POPULATION IS CONVERGED INTO A SUFFICIENTLY SMALL SPACE
IF (IPCNVG .EQ. 1) GO TO 9200
C
C NONE OF THE STOPPING CRITERIA IS SATISFIED, CONTINUE SEARCH
C
C CHECK FOR COMPLEX NUMBER REDUCTION
IF (NGS1 .GT .MINGS) THEN
NGS2 = NGS1
NGS1 = NGS1 - 1
NPT1 = NGS1 * NPG
CALL COMP(NOPT,NPT1,NGS1,NGS2,NPG,X,XF,CX,CF)
END IF
C
C END OF MAIN LOOP -----------
GO TO 1000
C
C SEARCH TERMINATED
9000 CONTINUE
WRITE(ISCE,800) MAXN,LOOP,IGS,NLOOP
WRITE(*,800) MAXN,LOOP,IGS,NLOOP
GO TO 9999
9100 CONTINUE
WRITE(ISCE,810) PCENTO*100.,KSTOP
WRITE(*,810) PCENTO*100.,KSTOP
GO TO 9999
9200 WRITE(ISCE,820) GNRNG*100.
WRITE(*,820) GNRNG*100.
9999 CONTINUE
C
C PRINT THE FINAL PARAMETER ESTIMATE AND ITS FUNCTION VALUE
WRITE(ISCE,830)
WRITE(ISCE,510) (XNAME(J),J=1,NOPT)
WRITE(ISCE,520) BESTF,(BESTX(J),J=1,NOPT)
WRITE(*,830)
WRITE(*,510) (XNAME(J),J=1,NOPT)
WRITE(*,520) BESTF,(BESTX(J),J=1,NOPT)
AF = BESTF
DO J = 1, NOPT
A(J) = BESTX(J)
END DO
10000 CONTINUE
C
C END OF SUBROUTINE SCEUA
RETURN
500 FORMAT(//,'*** PRINT THE INITIAL POINT AND ITS CRITERION ',
& 'VALUE ***')
510 FORMAT(/,' CRITERION',16(3X,A4,3X),/1X,80(1H-))
520 FORMAT(F8.3,17G10.3)
600 FORMAT(//,1X,'*** PRINT THE RESULTS OF THE SCE SEARCH ***')
610 FORMAT(/,1X,'LOOP',1X,'TRIALS',1X,'COMPLXS',1X,'BEST F',1X,
& 'WORST F',1X,'PAR RNG',1X,16(3X,A4,3X))
630 FORMAT(I5,1X,I5,3X,I5,3F8.3,17(G10.3))
650 FORMAT(/,1X,'POPULATION AT LOOP ',I3,/,1X,22(1H-))
660 FORMAT(F8.3,17(G10.3))
800 FORMAT(//,1X,'*** OPTIMIZATION SEARCH TERMINATED BECAUSE THE',
& ' LIMIT ON THE MAXIMUM',/,5X,'NUMBER OF TRIALS ',I5,
& ' EXCEEDED. SEARCH WAS STOPPED AT',/,5X,'SUB-COMPLEX ',
& I3,' OF COMPLEX ',I3,' IN SHUFFLING LOOP ',I3,' ***')
810 FORMAT(//,1X,'*** OPTIMIZATION TERMINATED BECAUSE THE CRITERION',
& ' VALUE HAS NOT CHANGED ',/,5X,F5.2,' PERCENT IN',I3,
& ' SHUFFLING LOOPS ***')
820 FORMAT(//,1X,'*** OPTIMIZATION TERMINATED BECAUSE THE POPULATION',
& ' HAS CONVERGED INTO ',/,4X,F5.2,' PERCENT OF THE',
& ' FEASIBLE SPACE ***')
830 FORMAT(//,'*** PRINT THE FINAL PARAMETER ESTIMATE AND ITS',
& ' CRITERION VALUE ***')
END
C
C
C
C====================================================================
SUBROUTINE CCE(NOPT,NPS,S,SF,BL,BU,XNSTD,ICALL,MAXN,ISEED)
C
C ALGORITHM GENERATE A NEW POINT(S) FROM A SUB-COMPLEX
C
C SUB-COMPLEX VARIABLES
PARAMETER (C1=0.8,C2=0.4)
DIMENSION S(50,16),SF(50),BU(16),BL(16),XNSTD(16)
C
C LIST OF LOCAL VARIABLES
C SB(.) = THE BEST POINT OF THE SIMPLEX
C SW(.) = THE WORST POINT OF THE SIMPLEX
C W2(.) = THE SECOND WORST POINT OF THE SIMPLEX
C FW = FUNCTION VALUE OF THE WORST POINT
C CE(.) = THE CENTROID OF THE SIMPLEX EXCLUDING WO
C SNEW(.) = NEW POINT GENERATED FROM THE SIMPLEX
C IVIOL = FLAG INDICATING IF CONSTRAINTS ARE VIOLATED
C = 1 , YES
C = 0 , NO
C
DIMENSION SW(16),SB(16),CE(16),SNEW(16)
C
C EQUIVALENCE OF VARIABLES FOR READABILTY OF CODE
N = NPS
M = NOPT
ALPHA = 1.0
BETA = 0.5
C
C IDENTIFY THE WORST POINT WO OF THE SUB-COMPLEX S
C COMPUTE THE CENTROID CE OF THE REMAINING POINTS
C COMPUTE STEP, THE VECTOR BETWEEN WO AND CE
C IDENTIFY THE WORST FUNCTION VALUE FW
DO J = 1, M
SB(J) = S(1,J)
SW(J) = S(N,J)
CE(J) = 0.0
DO I = 1, N-1
CE(J) = CE(J) + S(I,J)
END DO
CE(J) = CE(J)/DBLE(N-1)
END DO
FW = SF(N)
C
C COMPUTE THE NEW POINT SNEW
C
C FIRST TRY A REFLECTION STEP
DO J = 1, M
SNEW(J) = CE(J) + ALPHA * (CE(J) - SW(J))
END DO
C
C CHECK IF SNEW SATISFIES ALL CONSTRAINTS
CALL CHKCST(NOPT,SNEW,BL,BU,IBOUND)
C
C
C SNEW IS OUTSIDE THE BOUND,
C CHOOSE A POINT AT RANDOM WITHIN FEASIBLE REGION ACCORDING TO
C A NORMAL DISTRIBUTION WITH BEST POINT OF THE SUB-COMPLEX
C AS MEAN AND STANDARD DEVIATION OF THE POPULATION AS STD
IF (IBOUND .GE. 1) CALL GETPNT(NOPT,2,ISEED,SNEW,BL,BU,XNSTD,SB)
C
C
C COMPUTE THE FUNCTION VALUE AT SNEW
FNEW = FUNCTN(NOPT,SNEW)
ICALL = ICALL + 1
C
C COMPARE FNEW WITH THE WORST FUNCTION VALUE FW
C
C FNEW IS LESS THAN FW, ACCEPT THE NEW POINT SNEW AND RETURN
IF (FNEW .LE. FW) GO TO 2000
IF (ICALL .GE. MAXN) GO TO 3000
C
C
C FNEW IS GREATER THAN FW, SO TRY A CONTRACTION STEP
DO J = 1, M
SNEW(J) = CE(J) - BETA * (CE(J) - SW(J))
END DO
C
C COMPUTE THE FUNCTION VALUE OF THE CONTRACTED POINT
FNEW = FUNCTN(NOPT,SNEW)
ICALL = ICALL + 1
C
C COMPARE FNEW TO THE WORST VALUE FW
C IF FNEW IS LESS THAN OR EQUAL TO FW, THEN ACCEPT THE POINT AND RETURN
IF (FNEW .LE. FW) GO TO 2000
IF (ICALL .GE. MAXN) GO TO 3000
C
C
C IF BOTH REFLECTION AND CONTRACTION FAIL, CHOOSE ANOTHER POINT
C ACCORDING TO A NORMAL DISTRIBUTION WITH BEST POINT OF THE SUB-COMPLEX
C AS MEAN AND STANDARD DEVIATION OF THE POPULATION AS STD
1000 CALL GETPNT(NOPT,2,ISEED,SNEW,BL,BU,XNSTD,SB)
C
C COMPUTE THE FUNCTION VALUE AT THE RANDOM POINT
FNEW = FUNCTN(NOPT,SNEW)
ICALL = ICALL + 1
C
C
C REPLACE THE WORST POINT BY THE NEW POINT
2000 CONTINUE
DO J = 1, M
S(N,J) = SNEW(J)
END DO
SF(N) = FNEW
3000 CONTINUE
C
C END OF SUBROUTINE CCE
RETURN
END
C
C
C
C===================================================================
SUBROUTINE GETPNT(NOPT,IDIST,ISEED,X,BL,BU,STD,XI)
C
C THIS SUBROUTINE GENERATES A NEW POINT WITHIN FEASIBLE REGION
C
C X(.) = NEW POINT
C XI(.) = FOCAL POINT
C BL(.) = LOWER BOUND
C BU(.) = UPPER BOUND
C STD(.) = STANDARD DEVIATION OF PROBABILITY DISTRIBUTION
C IDIST = PROBABILITY FLAG
C = 1 - UNIFORM DISTRIBUTION
C = 2 - GAUSSIAN DISTRIBUTION
C
DIMENSION X(16),BL(16),BU(16),STD(16),XI(16)
C
1 DO J=1, NOPT
2 IF (IDIST .EQ. 1) RAND = RAN1(ISEED)
IF (IDIST .EQ. 2) RAND = GASDEV(ISEED)
X(J) = XI(J) + STD(J) * RAND * (BU(J) - BL(J))
C
C CHECK EXPLICIT CONSTRAINTS
C
CALL CHKCST(1,X(J),BL(J),BU(J),IBOUND)
IF (IBOUND .GE. 1) GO TO 2
END DO
C
C CHECK IMPLICIT CONSTRAINTS
C
CALL CHKCST(NOPT,X,BL,BU,IBOUND)
IF (IBOUND .GE. 1) GO TO 1
C
RETURN
END
C
C
C
C===================================================================
SUBROUTINE PARSTT(NPT,NOPT,X,XNSTD,BOUND,GNRNG,IPCNVG)
C
C SUBROUTINE CHECKING FOR PARAMETER CONVERGENCE
DIMENSION X(2000,16),XMAX(16),XMIN(16)
DIMENSION XMEAN(16),XNSTD(16),BOUND(16)
PARAMETER (DELTA = 1.0D-20,PEPS=1.0D-3)
C
C COMPUTE MAXIMUM, MINIMUM AND STANDARD DEVIATION OF PARAMETER VALUES
GSUM = 0.D0
DO K = 1, NOPT
XMAX(K) = -1.0D+20
XMIN(K) = 1.0D+20
XSUM1 = 0.D0
XSUM2 = 0.D0
DO I = 1, NPT
XMAX(K) = AMAX1(X(I,K), XMAX(K))
XMIN(K) = AMIN1(X(I,K), XMIN(K))
XSUM1 = XSUM1 + X(I,K)
XSUM2 = XSUM2 + X(I,K)*X(I,K)
END DO
XMEAN(K) = XSUM1 / NPT
XNSTD(K) = (XSUM2 / NPT - XMEAN(K)*XMEAN(K))
IF (XNSTD(K) .LE. DELTA) XNSTD(K) = DELTA
XNSTD(K) = SQRT(XNSTD(K))
XNSTD(K) = XNSTD(K) / BOUND(K)
GSUM = GSUM + LOG( DELTA + (XMAX(K)-XMIN(K))/BOUND(K) )
END DO
GNRNG = EXP(GSUM/NOPT)
C
C CHECK IF NORMALIZED STANDARD DEVIATION OF PARAMETER IS <= EPS
IPCNVG = 0
IF (GNRNG .LE. PEPS) THEN
IPCNVG = 1
END IF
C
C END OF SUBROUTINE PARSTT
RETURN
END
C
C
C
C====================================================================
SUBROUTINE COMP(N,NPT,NGS1,NGS2,NPG,A,AF,B,BF)
C
C
C THIS SUBROUTINE REDUCE INPUT MATRIX A(N,NGS2*NPG) TO MATRIX
C B(N,NGS1*NPG) AND VECTOR AF(NGS2*NPG) TO VECTOR BF(NGS1*NPG)
DIMENSION A(2000,16),AF(2000),B(2000,16),BF(2000)
DO IGS=1, NGS1
DO IPG=1, NPG
K1=(IPG-1)*NGS2 + IGS
K2=(IPG-1)*NGS1 + IGS
DO I=1, N
B(K2,I) = A(K1,I)
END DO
BF(K2) = AF(K1)
END DO
END DO
C
DO J=1, NPT
DO I=1, N
A(J,I) = B(J,I)
END DO
AF(J) = BF(J)
END DO
C
C END OF SUBROUTINE COMP
RETURN
END
C
C
C
C===================================================================
SUBROUTINE SORT(N,M,RB,RA)
C
C
C SORTING SUBROUTINE ADAPTED FROM "NUMERICAL RECIPES"
C BY W.H. PRESS ET AL., PP. 233-234
C
C LIST OF VARIABLES
C RA(.) = ARRAY TO BE SORTED
C RB(.,.) = ARRAYS ORDERED CORRESPONDING TO REARRANGEMENT OF RA(.)
C WK(.,.), IWK(.) = LOCAL VARIBLES
C
DIMENSION RA(2000),RB(2000,16),WK(2000,16),IWK(2000)
C
CALL INDEXX(N, RA, IWK)
DO 11 I = 1, N
WK(I,1) = RA(I)
11 CONTINUE
DO 12 I = 1, N
RA(I) = WK(IWK(I),1)
12 CONTINUE
DO 14 J = 1, M
DO 13 I = 1, N
WK(I,J) = RB(I,J)
13 CONTINUE
14 CONTINUE
DO 16 J = 1, M
DO 15 I = 1, N
RB(I,J) = WK(IWK(I),J)
15 CONTINUE
16 CONTINUE
C
C END OF SUBROUTINE SORT
RETURN
END
C
C
C===========================================================
SUBROUTINE SORT1(N,RA)
C
C
C SORTING SUBROUTINE ADAPTED FROM "NUMERICAL RECIPES"
C BY W.H. PRESS ET AL., PP. 231
C
C LIST OF VARIABLES
C RA(.) = INTEGER ARRAY TO BE SORTED
C
DIMENSION RA(N)
C
INTEGER RA, RRA
C
L = (N / 2) + 1
IR = N
10 CONTINUE
IF (L .GT. 1) THEN
L = L - 1
RRA = RA(L)
ELSE
RRA = RA(IR)
RA(IR) = RA(1)
IR = IR - 1
IF (IR .EQ. 1) THEN
RA(1) = RRA
RETURN
END IF
END IF
I = L
J = L + L
20 IF (J .LE. IR) THEN
IF (J .LT. IR) THEN
IF (RA(J) .LT. RA(J + 1)) J = J + 1
END IF
IF (RRA .LT. RA(J)) THEN
RA(I) = RA(J)
I = J
J = J + J
ELSE
J = IR + 1
END IF
GOTO 20
END IF
RA(I) = RRA
GOTO 10
C
C END OF SUBROUTINE SORT1
END
C
C
C
C=======================================================
SUBROUTINE INDEXX(N, ARRIN, INDX)
C
C
C THIS SUBROUTINE IS FROM "NUMERICAL RECIPES" BY PRESS ET AL.
DIMENSION ARRIN(N), INDX(N)
C
DO 11 J = 1, N
INDX(J) = J
11 CONTINUE
L = (N / 2) + 1
IR = N
10 CONTINUE
IF (L .GT. 1) THEN
L = L - 1
INDXT = INDX(L)
Q = ARRIN(INDXT)
ELSE
INDXT = INDX(IR)
Q = ARRIN(INDXT)
INDX(IR) = INDX(1)
IR = IR - 1
IF (IR .EQ. 1) THEN
INDX(1) = INDXT
RETURN
END IF
END IF
I = L
J = L + L
20 IF (J .LE. IR) THEN
IF (J .LT. IR) THEN
IF (ARRIN(INDX(J)) .LT. ARRIN(INDX(J + 1))) J = J + 1
END IF
IF (Q .LT. ARRIN(INDX(J))) THEN
INDX(I) = INDX(J)
I = J
J = J + J
ELSE
J = IR + 1
END IF
GOTO 20
END IF
INDX(I) = INDXT
GOTO 10
C
C END OF SUBROUTINE INDEXX
END
C
C
C
C==============================================================
REAL FUNCTION RAN1(IDUM)
C
C
C THIS SUBROUTINE IS FROM "NUMERICAL RECIPES" BY PRESS ET AL.
DIMENSION R(97)
PARAMETER (M1 = 259200, IA1 = 7141, IC1 = 54773, RM1 =
&3.8580247E-6)
PARAMETER (M2 = 134456, IA2 = 8121, IC2 = 28411, RM2 =
&7.4373773E-6)
PARAMETER (M3 = 243000, IA3 = 4561, IC3 = 51349)
SAVE
DATA IFF / 0 /
IF ((IDUM .LT. 0) .OR. (IFF .EQ. 0)) THEN
IFF = 1
IX1 = MOD(IC1 - IDUM,M1)
IX1 = MOD((IA1 * IX1) + IC1,M1)
IX2 = MOD(IX1,M2)
IX1 = MOD((IA1 * IX1) + IC1,M1)
IX3 = MOD(IX1,M3)
DO 11 J = 1, 97
IX1 = MOD((IA1 * IX1) + IC1,M1)
IX2 = MOD((IA2 * IX2) + IC2,M2)
R(J) = (DBLE(IX1) + (DBLE(IX2) * RM2)) * RM1
11 CONTINUE
IDUM = 1
END IF
IX1 = MOD((IA1 * IX1) + IC1,M1)
IX2 = MOD((IA2 * IX2) + IC2,M2)
IX3 = MOD((IA3 * IX3) + IC3,M3)
J = 1 + ((97 * IX3) / M3)
IF ((J .GT. 97) .OR. (J .LT. 1)) PAUSE
RAN1 = R(J)
R(J) = (DBLE(IX1) + (DBLE(IX2) * RM2)) * RM1
C
C END OF SUBROUTINE RAN1
RETURN
END
C
C
C
C===============================================================
FUNCTION GASDEV(IDUM)
C
C
C THIS SUBROUTINE IS FROM "NUMERICAL RECIPES" BY PRESS ET AL.
COMMON /GASBLK/ ISET
DATA ISET / 0 /
IF (ISET .EQ. 0) THEN
1 V1 = (2. * RAN1(IDUM)) - 1.
V2 = (2. * RAN1(IDUM)) - 1.
R = (V1 ** 2) + (V2 ** 2)
IF (R .GE. 1.) GOTO 1
FAC = SQRT(- ((2. * LOG(R)) / R))
GSET = V1 * FAC
GASDEV = V2 * FAC
ISET = 1
ELSE
GASDEV = GSET
ISET = 0
END IF
C
C END OF SUBROUTINE GASDEV
RETURN
END
scein.for：
SUBROUTINE SCEIN(NOPT,MAXN,KSTOP,PCENTO,ISEED,
& NGS,NPG,NPS,NSPL,MINGS,INIFLG,IPRINT)
C
C
C THIS SUBROUTINE READS AND PRINTS THE INPUT VARIABLES FOR
C SHUFFLED COMPLEX EVOLUTION METHOD FOR GLOBAL OPTIMIZATION
C -- VERSION 2.1
C
C WRITTEN BY QINGYUN DUAN - UNIVERSITY OF ARIZONA, APRIL 1992
C
C
CHARACTER*10 PCNTRL,DEFLT,USRSP
CHARACTER*4 REDUC,INITL,YSFLG,NOFLG
COMMON/IOPARS/ICNTRL,IOUT,IDAT,IWBAL,ISCE,IPE,IPC,IDET
C
DATA DEFLT/' DEFAULT '/
DATA USRSP/'USER SPEC.'/
DATA YSFLG/'YES '/
DATA NOFLG/'NO '/
C
IERROR = 0
IWARN = 0
WRITE(ISCE,700)
700 FORMAT(10X,'SHUFFLED COMPLEX EVOLUTION GLOBAL OPTIMIZATION',
& /,10X,46(1H=))
C
C
C READ THE SCE CONTROL PARAMETERS
IDEFLT = 0
READ (ICNTRL,*)
READ (ICNTRL,800) MAXN,KSTOP,PCENTO,NGS,ISEED,IDEFLT
800 FORMAT(2I5,F5.2,3I5)
IF (ISEED .EQ. 0) ISEED = 1969
C
C
C IF IDEFLT IS EQUAL TO 1, READ THE SCE CONTROL PARAMETERS
IF (IDEFLT .EQ. 1) THEN
READ(ICNTRL,810) NPG,NPS,NSPL,MINGS,INIFLG,IPRINT
810 FORMAT(6I5)
PCNTRL = USRSP
IF (NPG .EQ. 0) NPG = 2*NOPT + 1
IF (NPS .EQ. 0) NPS = NOPT + 1
IF (NSPL .EQ. 0) NSPL = NPG
IF (MINGS .EQ. 0) MINGS = NGS
ELSE
READ(ICNTRL,*)
PCNTRL = DEFLT
END IF
C
C
C IF IDEFLT IS EQUAL TO 0, SET THE SCE CONTROL PARAMETERS TO
C THE DEFAULT VALUES
IF (IDEFLT .EQ. 0) THEN
NPG = 2*NOPT + 1
NPS = NOPT + 1
NSPL = NPG
MINGS = NGS
INIFLG = 0
IPRINT = 0
END IF
C
C
C CHECK IF THE SCE CONTROL PARAMETERS ARE VALID
IF (NGS .LT. 1 .OR. NGS .GE. 1320) THEN
WRITE(ISCE,900) NGS
900 FORMAT(//,1X,'**ERROR** NUMBER OF COMPLEXES IN INITIAL ',
* ' POPULATION ',I5,' IS NOT A VALID CHOICE')
IERROR = IERROR + 1
END IF
C
IF (KSTOP .LT. 0 .OR. KSTOP .GE. 10) THEN
WRITE(ISCE,901) KSTOP
901 FORMAT(//,1X,'**WARNING** THE NUMBER OF SHUFFLING LOOPS IN',
* ' WHICH THE CRITERION VALUE MUST CHANGE ',/,13X,'SHOULD BE',
* ' GREATER THAN 0 AND LESS THAN 10. ','KSTOP = ',I2,
* ' WAS SPECIFIED.'/,13X,'BUT KSTOP = 5 WILL BE USED INSTEAD.')
IWARN = IWARN + 1
KSTOP=5
END IF
C
IF (MINGS .LT. 1 .OR. MINGS .GT. NGS) THEN
WRITE(ISCE,902) MINGS
902 FORMAT(//,1X,'**WARNING** THE MINIMUM NUMBER OF COMPLEXES ',
* I2,' IS NOT A VALID CHOICE. SET IT TO DEFAULT')
IWARN = IWARN + 1
MINGS = NGS
END IF
C
IF (NPG .LT. 2 .OR. NPG .GT. 1320/MAX(NGS,1)) THEN
WRITE(ISCE,903) NPG
903 FORMAT(//,1X,'**WARNING** THE NUMBER OF POINTS IN A COMPLEX ',
* I4,' IS NOT A VALID CHOICE, SET IT TO DEFAULT')
IWARN = IWARN + 1
NPG = 2*NOPT+1
END IF
C
IF (NPS.LT.2 .OR. NPS.GT.NPG .OR. NPS.GT.50) THEN
WRITE(ISCE,904) NPS
904 FORMAT(//,1X,'**WARNING** THE NUMBER OF POINTS IN A SUB-',
* 'COMPLEX ',I4,' IS NOT A VALID CHOICE, SET IT TO DEFAULT')
IWARN = IWARN + 1
NPS = NOPT + 1
END IF
C
IF (NSPL .LT. 1) THEN
WRITE(ISCE,905) NSPL
905 FORMAT(//,1X,'**WARNING** THE NUMBER OF EVOLUTION STEPS ',
* 'TAKEN IN EACH COMPLEX BEFORE SHUFFLING ',I4,/,13X,
* 'IS NOT A VALID CHOICE, SET IT TO DEFAULT')
IWARN = IWARN + 1
NSPL = NPG
END IF
C
C COMPUTE THE TOTAL NUMBER OF POINTS IN INITIAL POPULATION
NPT = NGS * NPG
C
IF (NPT .GT. 1320) THEN
WRITE(ISCE,906) NPT
906 FORMAT(//,1X,'**WARNING** THE NUMBER OF POINTS IN INITIAL ',
* 'POPULATION ',I5,' EXCEED THE POPULATION LIMIT,',/,13X,
* 'SET NGS TO 2, AND NPG, NPS AND NSPL TO DEFAULTS')
IWARN = IWARN + 1
NGS = 2
NPG = 2*NOPT + 1
NPS = NOPT + 1
NSPL = NPG
END IF
C
C PRINT OUT THE TOTAL NUMBER OF ERROR AND WARNING MESSAGES
IF (IERROR .GE. 1) WRITE(ISCE,907) IERROR
907 FORMAT(//,1X,'*** TOTAL NUMBER OF ERROR MESSAGES IS ',I2)
C
IF (IWARN .GE. 1) WRITE(ISCE,908) IWARN
908 FORMAT(//,1X,'*** TOTAL NUMBER OF WARNING MESSAGES IS ',I2)
C
IF (MINGS .LT. NGS) THEN
REDUC = YSFLG
ELSE
REDUC = NOFLG
END IF
C
IF (INIFLG .NE. 0) THEN
INITL = YSFLG
ELSE
INITL = NOFLG
END IF
C
C
C PRINT SHUFFLED COMPLEX EVOLUTION OPTIMIZATION OPTIONS
104 WRITE(ISCE,910)
910 FORMAT(//,2X,'SCE CONTROL',5X,'MAX TRIALS',5X,
&'REQUIRED IMPROVEMENT',5X,'RANDOM',/,3X,'PARAMETER',8X,
&'ALLOWED',6X,'PERCENT',4X,'NO. LOOPS',6X,'SEED',/,
&2X,11(1H-),5X,10(1H-),5X,7(1H-),4X,9(1H-),5X,6(1H-))
C
PCENTA=PCENTO*100.
WRITE(ISCE,912) PCNTRL,MAXN,PCENTA,KSTOP,ISEED
912 FORMAT(3X,A10,7X,I5,10X,F3.1,9X,I2,9X,I5)
WRITE(ISCE,914) NGS,NPG,NPT,NPS,NSPL
914 FORMAT(//,18X,'SCE ALGORITHM CONTROL PARAMETERS',/,18X,32(1H=),
&//,2X,'NUMBER OF',5X,'POINTS PER',5X,'POINTS IN',6X,'POINTS PER',
&4X,'EVOL. STEPS',/,2X,'COMPLEXES',6X,'COMPLEX',6X,'INI. POPUL.',
&5X,'SUB-COMPLX',4X,'PER COMPLEX',/,2X,9(1H-),5X,10(1H-),4X,
&11(1H-),5X,10(1H-),4X,11(1H-),5X,/,2X,5(I5,10X))
WRITE(ISCE,915) REDUC,MINGS,INITL
915 FORMAT(//,15X,'COMPLX NO.',5X,'MIN COMPLEX',5X,'INI. POINT',/,
&15X,'REDUCTION',6X,'NO. ALLOWED',6X,'INCLUDED',/,
&15X,10(1H-),5X,11(1H-),5X,10(1H-),/,18X,A4,6X,I8,13X,A4)
IF (IERROR .GE. 1) THEN
WRITE(ISCE,922)
922 FORMAT(//,'*** THE OPTIMIZATION SEARCH IS NOT CONDUCTED BECAUSE',
& ' OF INPUT DATA ERROR ***')
STOP
END IF
C
C END OF SUBROUTINE SCEIN
RETURN
END
翻译
以下代码仅供参考，都无法成功运行。
笔者根据源码翻译如下：
sceua.for
SUBROUTINE SCEUA(A,AF,BL,BU,NOPT,MAXN,KSTOP,PCENTO,ISEED,
& NGS,NPG,NPS,NSPL,MINGS,INIFLG,IPRINT)
C
C
C 洗牌复形演化算法（全局优化）
C ——版本2.1
C
C 作者：段青云
C 亚利桑那州图森市亚利桑那大学
C 水文与水资源系 AZ 85721
C (602) 621-9360, 邮箱: DUAN@HWR.ARIZONA.EDU
C
C 写于1990年10月
C 1991年8月修改
C 1992年4月修改
C
C 作者声明：
C --------------------
C
C 该通用全局优化程序由亚利桑那大学水文与水资源系开发。
C 有关SCE-UA方法的更多信息
C 可以根据上面列出的地址和电话号码
C 咨询段青云博士、S. SOROOSHIAN博士或V.K. GUPTA博士。
C 我们要求本程序的所有用户正确引用论文：
C Duan, Q., Sorooshian, S., and Gupta, V. (1992),
C Effective and efficient global optimization
C for conceptual rainfall-runoff models,
C Water Resour. Res., 28(4), 1015–1031, doi:10.1029/91WR02985.
C
C
C 输入参数变量列表
C
C A(.) = 初始参数集
C BL(.) = 参数下限
C BU(.) = 参数上限
C NOPT = 待优化的参数数量
C
C
C SCE算法控制参数列表：
C
C NGS = 初始种群中的复形数量
C NPG = 每个复形中的点的个数
C NPT = 初始种群中的总点数(NPT=NGS*NPG)
C NPS = 子复形中的点数
C NSPL = 复形洗牌前每个复形
C 允许的进化步骤数
C MINGS = 如果在优化过程中允许减少复形的数量，
C 则需要最少的复形数量
C ISEED = 初始随机种子
C INIFLG = 标记是否将初始点包括在种群中
C = 0, 不包括
C = 1, 包括
C IPRINT = 用于控制每个洗牌循环后的打印输出的标志
C = 0, 打印关于种群最佳点的信息
C = 1, 打印关于种群每个点的信息
C
C
C 收敛检查参数
C
C MAXN = 优化终止前允许的最大试验次数
C KSTOP = 在终止优化之前，
C 标准值必须改变给定百分比的洗牌循环数
C PCENTO = 给定数量的随机循环中，
C 标准值必须改变的百分比
C IPCNVG = 指示是否达到参数收敛的标志
C （即，检查 GNRNG 是否小于 0.001）
C = 0, 参数收敛不满足
C = 1, 参数收敛满足
C
C
C 局部变量列表
C X(.,.) = 种群中点的坐标
C XF(.) = X(.,.)的函数值
C XX(.) = X中单点的坐标
C CX(.,.) = 复形中点的坐标
C CF(.) = CX(.,.)的函数值
C S(.,.) = 当前单纯形中点的坐标
C SF(.) = S(.,.)的函数值
C BESTX(.) = 当前洗牌循环的最佳点
C BESTF = BESTX(.)的函数值
C WORSTX(.) = 当前洗牌循环的最坏点
C WORSTF = WORSTX(.)的函数值
C XNSTD(.) = 种群中参数的标准差
C GNRNG = 参数范围的归一化几何平均值
C LCS(.) = 索引定位S(.,.)在 X(.,.)中的位置
C BOUND(.) = 优化变量的约束
C NGS1 = 当前种群中的复形数量
C NGS2 = 前一种群中的复形数量
C ISEED1 = 当前随机种子
C CRITER(.) = 包含最近10个洗牌循环的
C 最佳标准值的向量
C
CHARACTER*4 XNAME(16)
C
C 来自输入数据的数组
DIMENSION A(16),BL(16),BU(16)
C
C 局部数组
DIMENSION X(2000,16),XX(16),BESTX(16),WORSTX(16),XF(2000)
DIMENSION S(50,16),SF(50),LCS(50),CX(2000,16),CF(2000)
DIMENSION XNSTD(16),BOUND(16),CRITER(20),UNIT(16)
C
COMMON/IOPARS/ICNTRL,IOUT,IDAT,IWBAL,ISCE,IPE,IPC,IDET
C
DATA XNAME /' X1',' X2',' X3',' X4',' X5',' X6',' X7',
&' X8',' X9',' X10',' X11',' X12',' X13',' X14',' X15',' X16'/
C
C 初始化变量
WRITE(ISCE,400)
400 FORMAT(//,2X,50(1H=),/,2X,'ENTER THE SHUFFLED COMPLEX EVOLUTION',
& ' GLOBAL SEARCH',/,2X,50(1H=))
WRITE(*,400)
NLOOP = 0
LOOP = 0
IGS = 0
C
C 初始化随机种子为负整数
ISEED1 = -ABS(ISEED)
C
C 计算在初始种群中的总点数
NPT = NGS * NPG
NGS1 = NGS
NPT1 = NPT
C
C 计算待优化参数约束范围
DO J = 1, NOPT
BOUND(J) = BU(J) - BL(J)
UNIT(J) = 1.0
END DO
C
C 计算初始点函数值
FA = FUNCTN(NOPT,A)
C
C 打印初始点及其标准值
WRITE(ISCE,500)
WRITE(*, 500)
WRITE(ISCE,510) (XNAME(J),J=1,NOPT)
WRITE(*, 510) (XNAME(J),J=1,NOPT)
WRITE(ISCE,520) FA,(A(J),J=1,NOPT)
WRITE(*, 520) FA,(A(J),J=1,NOPT)
IF (MAXN .EQ. 1) GO TO 10000
C
C 在NPT1点的参数空间内生成一个初始集合
C 如果INIFLG等于1，设置X(1,.)为初始点A(.)
IF (INIFLG .EQ. 1) THEN
DO J = 1, NOPT
X(1,J) = A(J)
END DO
XF(1) = FA
C
C 否则，随机生成一个点并设置其等于X(1,.)
ELSE
CALL GETPNT(NOPT,1,ISEED1,XX,BL,BU,UNIT,BL)
DO J=1, NOPT
X(1,J) = XX(J)
END DO
XF(1) = FUNCTN(NOPT,XX)
END IF
ICALL = 1
IF (ICALL .GE. MAXN) GO TO 9000
C
C 生成在参数空间均匀分布的NPT1-1个随机点，
C 并计算对应的函数值
DO I = 2, NPT1
CALL GETPNT(NOPT,1,ISEED1,XX,BL,BU,UNIT,BL)
DO J = 1, NOPT
X(I,J) = XX(J)
END DO
XF(I) = FUNCTN(NOPT,XX)
ICALL = ICALL + 1
IF (ICALL .GE. MAXN) THEN
NPT1 = I
GO TO 45
END IF
END DO
C
C 按函数值递增的顺序排列点
45 CALL SORT(NPT1,NOPT,X,XF)
C
C 记录最好和最差的点
DO J = 1, NOPT
BESTX(J) = X(1,J)
WORSTX(J) = X(NPT1,J)
END DO
BESTF = XF(1)
WORSTF = XF(NPT1)
C
C 计算初始种群的参数范围
CALL PARSTT(NPT1,NOPT,X,XNSTD,BOUND,GNRNG,IPCNVG)
C
C 打印初始种群的结果
WRITE(ISCE,600)
WRITE(*, 600)
WRITE(ISCE,610) (XNAME(J),J=1,NOPT)
WRITE(*, 610) (XNAME(J),J=1,NOPT)
WRITE(ISCE,630) NLOOP,ICALL,NGS1,BESTF,WORSTF,GNRNG,
& (BESTX(J),J=1,NOPT)
WRITE(*, 630) NLOOP,ICALL,NGS1,BESTF,WORSTF,GNRNG,
& (BESTX(J),J=1,NOPT)
IF (IPRINT .EQ. 1) THEN
WRITE(ISCE,650) NLOOP
DO I = 1, NPT1
WRITE(ISCE,660) XF(I),GNRNG,(X(I,J),J=1,NOPT)
END DO
END IF
C
IF (ICALL .GE. MAXN) GO TO 9000
IF (IPCNVG .EQ. 1) GO TO 9200
C
C 开始主循环 ----------------
1000 CONTINUE
NLOOP = NLOOP + 1
C
C 在复形上开始循环
DO 3000 IGS = 1, NGS1
C
C 为复形分配点
DO K1 = 1, NPG
K2 = (K1-1) * NGS1 + IGS
DO J = 1, NOPT
CX(K1,J) = X(K2,J)
END DO
CF(K1) = XF(K2)
END DO
C
C 开始内循环 - 子复形的随机选择 ---------------
DO 2000 LOOP = 1, NSPL
C
C 根据线性概率分布
C 选择子复形（NPS个点）
IF (NPS .EQ. NPG) THEN
DO K = 1, NPS
LCS(K) = K
END DO
GO TO 85
END IF
C
RAND = RAN1(ISEED1)
LCS(1) = 1 + INT(NPG + 0.5 - SQRT( (NPG+.5)**2 -
& NPG * (NPG+1) * RAND ))
DO K = 2, NPS
60 RAND = RAN1(ISEED1)
LPOS = 1 + INT(NPG + 0.5 - SQRT((NPG+.5)**2 -
& NPG * (NPG+1) * RAND ))
DO K1 = 1, K-1
IF (LPOS .EQ. LCS(K1)) GO TO 60
END DO
LCS(K) = LPOS
END DO
C
C 按照函数值递增的顺序排列子复形
CALL SORT1(NPS,LCS)
C
C 创建子复形数组
85 DO K = 1, NPS
DO J = 1, NOPT
S(K,J) = CX(LCS(K),J)
END DO
SF(K) = CF(LCS(K))
END DO
C
C 使用子复形生成新点
CALL CCE(NOPT,NPS,S,SF,BL,BU,XNSTD,ICALL,MAXN,ISEED1)
C
C 如果子复形被接受，
C 将新的子复形更换到复形中
DO K = 1, NPS
DO J = 1, NOPT
CX(LCS(K),J) = S(K,J)
END DO
CF(LCS(K)) = SF(K)
END DO
C
C 排序点
CALL SORT(NPG,NOPT,CX,CF)
C
C 如果超过最大运行次数，则跳出循环
IF (ICALL .GE. MAXN) GO TO 2222
C
C 结束内循环 ------------
2000 CONTINUE
2222 CONTINUE
C
C 将新复形替换为原始数组 X(.,.)
DO K1 = 1, NPG
K2 = (K1-1) * NGS1 + IGS
DO J = 1, NOPT
X(K2,J) = CX(K1,J)
END DO
XF(K2) = CF(K1)
END DO
IF (ICALL .GE. MAXN) GO TO 3333
C
C 结束复形循环
3000 CONTINUE
C
C 对点重新排序
3333 CALL SORT(NPT1,NOPT,X,XF)
C
C 记录最好和最坏的点
DO J = 1, NOPT
BESTX(J) = X(1,J)
WORSTX(J) = X(NPT1,J)
END DO
BESTF = XF(1)
WORSTF = XF(NPT1)
C
C 测试种群的参数收敛性
CALL PARSTT(NPT1,NOPT,X,XNSTD,BOUND,GNRNG,IPCNVG)
C
C 打印当前种群的结果
WRITE(ISCE,630) NLOOP,ICALL,NGS1,BESTF,WORSTF,GNRNG,
& (BESTX(J),J=1,NOPT)
WRITE(*,630) NLOOP,ICALL,NGS1,BESTF,WORSTF,GNRNG,
& (BESTX(J),J=1,NOPT)
IF (IPRINT .EQ. 1) THEN
WRITE(ISCE,650) NLOOP
DO I = 1, NPT1
WRITE(ISCE,660) XF(I),GNRNG,(X(I,J),J=1,NOPT)
END DO
END IF
C
C 测试是否超过了函数评估次数的最大值
IF (ICALL .GE. MAXN) GO TO 9000
C
C 计算无函数改进的连续循环次数
CRITER(20) = BESTF
IF (NLOOP .LT. (KSTOP+1)) GO TO 132
DENOMI = ABS(CRITER(20-KSTOP) + CRITER(20)) / 2.
TIMEOU = ABS(CRITER(20-KSTOP) - CRITER(20)) / DENOMI
IF (TIMEOU .LT. PCENTO) GO TO 9100
132 CONTINUE
DO L = 1, 19
CRITER(L) = CRITER(L+1)
END DO
C
C 如果种群聚集到一个足够小的空间
IF (IPCNVG .EQ. 1) GO TO 9200
C
C 不满足任何停止标准，继续搜索
C
C 检查复形数量减少
IF (NGS1 .GT .MINGS) THEN
NGS2 = NGS1
NGS1 = NGS1 - 1
NPT1 = NGS1 * NPG
CALL COMP(NOPT,NPT1,NGS1,NGS2,NPG,X,XF,CX,CF)
END IF
C
C 结束主循环 -----------
GO TO 1000
C
C 搜索终止
9000 CONTINUE
WRITE(ISCE,800) MAXN,LOOP,IGS,NLOOP
WRITE(*,800) MAXN,LOOP,IGS,NLOOP
GO TO 9999
9100 CONTINUE
WRITE(ISCE,810) PCENTO*100.,KSTOP
WRITE(*,810) PCENTO*100.,KSTOP
GO TO 9999
9200 WRITE(ISCE,820) GNRNG*100.
WRITE(*,820) GNRNG*100.
9999 CONTINUE
C
C 打印最终参数估计及其函数值
WRITE(ISCE,830)
WRITE(ISCE,510) (XNAME(J),J=1,NOPT)
WRITE(ISCE,520) BESTF,(BESTX(J),J=1,NOPT)
WRITE(*,830)
WRITE(*,510) (XNAME(J),J=1,NOPT)
WRITE(*,520) BESTF,(BESTX(J),J=1,NOPT)
AF = BESTF
DO J = 1, NOPT
A(J) = BESTX(J)
END DO
10000 CONTINUE
C
C 子程序SCEUA结束
RETURN
500 FORMAT(//,'*** PRINT THE INITIAL POINT AND ITS CRITERION ',
& 'VALUE ***')
510 FORMAT(/,' CRITERION',16(3X,A4,3X),/1X,80(1H-))
520 FORMAT(F8.3,17G10.3)
600 FORMAT(//,1X,'*** PRINT THE RESULTS OF THE SCE SEARCH ***')
610 FORMAT(/,1X,'LOOP',1X,'TRIALS',1X,'COMPLXS',1X,'BEST F',1X,
& 'WORST F',1X,'PAR RNG',1X,16(3X,A4,3X))
630 FORMAT(I5,1X,I5,3X,I5,3F8.3,17(G10.3))
650 FORMAT(/,1X,'POPULATION AT LOOP ',I3,/,1X,22(1H-))
660 FORMAT(F8.3,17(G10.3))
800 FORMAT(//,1X,'*** OPTIMIZATION SEARCH TERMINATED BECAUSE THE',
& ' LIMIT ON THE MAXIMUM',/,5X,'NUMBER OF TRIALS ',I5,
& ' EXCEEDED. SEARCH WAS STOPPED AT',/,5X,'SUB-COMPLEX ',
& I3,' OF COMPLEX ',I3,' IN SHUFFLING LOOP ',I3,' ***')
810 FORMAT(//,1X,'*** OPTIMIZATION TERMINATED BECAUSE THE CRITERION',
& ' VALUE HAS NOT CHANGED ',/,5X,F5.2,' PERCENT IN',I3,
& ' SHUFFLING LOOPS ***')
820 FORMAT(//,1X,'*** OPTIMIZATION TERMINATED BECAUSE THE POPULATION',
& ' HAS CONVERGED INTO ',/,4X,F5.2,' PERCENT OF THE',
& ' FEASIBLE SPACE ***')
830 FORMAT(//,'*** PRINT THE FINAL PARAMETER ESTIMATE AND ITS',
& ' CRITERION VALUE ***')
END
C
C
C
C====================================================================
SUBROUTINE CCE(NOPT,NPS,S,SF,BL,BU,XNSTD,ICALL,MAXN,ISEED)
C
C 该算法从子复形生成新点
C
C 子复形变量
PARAMETER (C1=0.8,C2=0.4)
DIMENSION S(50,16),SF(50),BU(16),BL(16),XNSTD(16)
C
C 局部变量列表
C SB(.) = 单纯形的最佳点
C SW(.) = 单纯形的最坏点WO（WORST POINT，简称WO）
C W2(.) = 单纯形的第二最坏点
C FW = 最坏点的函数值
C CE(.) = 单纯形排除最坏点（WO）的质心
C SNEW(.) = 从单纯形生成的新点
C IVIOL = 指示是否违反约束的标志
C = 1 , YES
C = 0 , NO
C
DIMENSION SW(16),SB(16),CE(16),SNEW(16)
C
C 代码可读性变量的等价性
N = NPS
M = NOPT
ALPHA = 1.0
BETA = 0.5
C
C 确定子复形S的最坏点WO
C 计算剩余点的质心CE（CENTROID，简称CE）
C 计算步长，最坏点WO和质心CE之间的向量
C 确定最坏的函数值FW（The worst function value，简称FW）
DO J = 1, M
SB(J) = S(1,J)
SW(J) = S(N,J)
CE(J) = 0.0
DO I = 1, N-1
CE(J) = CE(J) + S(I,J)
END DO
CE(J) = CE(J)/DBLE(N-1)
END DO
FW = SF(N)
C
C 计算新点SNEW
C
C 首先尝试反射步骤
DO J = 1, M
SNEW(J) = CE(J) + ALPHA * (CE(J) - SW(J))
END DO
C
C 检查是否SNEW满足所有约束
CALL CHKCST(NOPT,SNEW,BL,BU,IBOUND)
C
C
C SNEW在界外，
C 按正态分布在可行区域内随机选择一个点，
C 子复形的最佳点作为种群的均值
C 标准差作为STD
IF (IBOUND .GE. 1) CALL GETPNT(NOPT,2,ISEED,SNEW,BL,BU,XNSTD,SB)
C
C
C 计算SNEW点的函数值FNEW
FNEW = FUNCTN(NOPT,SNEW)
ICALL = ICALL + 1
C
C 将FNEW与最差函数值FW进行比较
C
C FNEW小于FW，接受新点SNEW并返回
IF (FNEW .LE. FW) GO TO 2000
IF (ICALL .GE. MAXN) GO TO 3000
C
C
C FNEW大于FW，所以尝试收缩步骤
DO J = 1, M
SNEW(J) = CE(J) - BETA * (CE(J) - SW(J))
END DO
C
C 计算收缩点的函数值FNEW
FNEW = FUNCTN(NOPT,SNEW)
ICALL = ICALL + 1
C
C 将FNEW与最差值FW进行比较
C 如果FNEW小于或等于FW，则接受该点并返回
IF (FNEW .LE. FW) GO TO 2000
IF (ICALL .GE. MAXN) GO TO 3000
C
C
C 如果反射和收缩都失败，
C 则根据正态分布选择另一个点，
C 子复形的最佳点作为种群的平均值并且标准偏差作为STD
1000 CALL GETPNT(NOPT,2,ISEED,SNEW,BL,BU,XNSTD,SB)
C
C 计算随机点的函数值
FNEW = FUNCTN(NOPT,SNEW)
ICALL = ICALL + 1
C
C
C 用新点替换最坏点
2000 CONTINUE
DO J = 1, M
S(N,J) = SNEW(J)
END DO
SF(N) = FNEW
3000 CONTINUE
C
C 子程序CCE结束
RETURN
END
C
C
C
C===================================================================
SUBROUTINE GETPNT(NOPT,IDIST,ISEED,X,BL,BU,STD,XI)
C
C 该子程序在可行区域内生成一个新点
C
C X(.) = 新点
C XI(.) = 焦点
C BL(.) = 下限
C BU(.) = 上限
C STD(.) = 概率分布的标准差
C IDIST = 概率标志
C = 1 - 均匀分布
C = 2 - 高斯分布即正态分布
C
DIMENSION X(16),BL(16),BU(16),STD(16),XI(16)
C
1 DO J=1, NOPT
2 IF (IDIST .EQ. 1) RAND = RAN1(ISEED)
IF (IDIST .EQ. 2) RAND = GASDEV(ISEED)
X(J) = XI(J) + STD(J) * RAND * (BU(J) - BL(J))
C
C 检查显式约束
C
CALL CHKCST(1,X(J),BL(J),BU(J),IBOUND)
IF (IBOUND .GE. 1) GO TO 2
END DO
C
C 检查隐式约束
C
CALL CHKCST(NOPT,X,BL,BU,IBOUND)
IF (IBOUND .GE. 1) GO TO 1
C
RETURN
END
C
C
C
C===================================================================
SUBROUTINE PARSTT(NPT,NOPT,X,XNSTD,BOUND,GNRNG,IPCNVG)
C
C 参数收敛的检查子程序
DIMENSION X(2000,16),XMAX(16),XMIN(16)
DIMENSION XMEAN(16),XNSTD(16),BOUND(16)
PARAMETER (DELTA = 1.0D-20,PEPS=1.0D-3)
C
C 计算参数值的最大值、最小值和标准差
GSUM = 0.D0
DO K = 1, NOPT
XMAX(K) = -1.0D+20
XMIN(K) = 1.0D+20
XSUM1 = 0.D0
XSUM2 = 0.D0
DO I = 1, NPT
XMAX(K) = AMAX1(X(I,K), XMAX(K))
XMIN(K) = AMIN1(X(I,K), XMIN(K))
XSUM1 = XSUM1 + X(I,K)
XSUM2 = XSUM2 + X(I,K)*X(I,K)
END DO
XMEAN(K) = XSUM1 / NPT
XNSTD(K) = (XSUM2 / NPT - XMEAN(K)*XMEAN(K))
IF (XNSTD(K) .LE. DELTA) XNSTD(K) = DELTA
XNSTD(K) = SQRT(XNSTD(K))
XNSTD(K) = XNSTD(K) / BOUND(K)
GSUM = GSUM + LOG( DELTA + (XMAX(K)-XMIN(K))/BOUND(K) )
END DO
GNRNG = EXP(GSUM/NOPT)
C
C 检查参数的标准化标准差是否 <= EPS
IPCNVG = 0
IF (GNRNG .LE. PEPS) THEN
IPCNVG = 1
END IF
C
C 子程序PARSTT结束
RETURN
END
C
C
C
C====================================================================
SUBROUTINE COMP(N,NPT,NGS1,NGS2,NPG,A,AF,B,BF)
C
C
C 此子程序将输入矩阵 A(N,NGS2*NPG) 减少到矩阵 B(N,NGS1*NPG)
C 并将向量 AF(NGS2*NPG) 减少到向量 BF(NGS1*NPG)
DIMENSION A(2000,16),AF(2000),B(2000,16),BF(2000)
DO IGS=1, NGS1
DO IPG=1, NPG
K1=(IPG-1)*NGS2 + IGS
K2=(IPG-1)*NGS1 + IGS
DO I=1, N
B(K2,I) = A(K1,I)
END DO
BF(K2) = AF(K1)
END DO
END DO
C
DO J=1, NPT
DO I=1, N
A(J,I) = B(J,I)
END DO
AF(J) = BF(J)
END DO
C
C 子程序COMP结束
RETURN
END
C
C
C
C===================================================================
SUBROUTINE SORT(N,M,RB,RA)
C
C
C 改编自“数字食谱”（"NUMERICAL RECIPES"）的排序子程序
C 由 W.H. PRESS等人，PP. 233-234
C
C 变量列表
C RA(.) = 要排序的数组
C RB(.,.) = 对应于 RA(.) 重排的阵列排序
C WK(.,.), IWK(.) = 局部变量
C
DIMENSION RA(2000),RB(2000,16),WK(2000,16),IWK(2000)
C
CALL INDEXX(N, RA, IWK)
DO 11 I = 1, N
WK(I,1) = RA(I)
11 CONTINUE
DO 12 I = 1, N
RA(I) = WK(IWK(I),1)
12 CONTINUE
DO 14 J = 1, M
DO 13 I = 1, N
WK(I,J) = RB(I,J)
13 CONTINUE
14 CONTINUE
DO 16 J = 1, M
DO 15 I = 1, N
RB(I,J) = WK(IWK(I),J)
15 CONTINUE
16 CONTINUE
C
C 子程序SORT结束
RETURN
END
C
C
C===========================================================
SUBROUTINE SORT1(N,RA)
C
C
C 改编自“数字食谱”（"NUMERICAL RECIPES"）的排序子程序
C 由 W.H. PRESS等人，PP. 231
C
C 变量列表
C RA(.) = 要排序的整数数组
C
DIMENSION RA(N)
C
INTEGER RA, RRA
C
L = (N / 2) + 1
IR = N
10 CONTINUE
IF (L .GT. 1) THEN
L = L - 1
RRA = RA(L)
ELSE
RRA = RA(IR)
RA(IR) = RA(1)
IR = IR - 1
IF (IR .EQ. 1) THEN
RA(1) = RRA
RETURN
END IF
END IF
I = L
J = L + L
20 IF (J .LE. IR) THEN
IF (J .LT. IR) THEN
IF (RA(J) .LT. RA(J + 1)) J = J + 1
END IF
IF (RRA .LT. RA(J)) THEN
RA(I) = RA(J)
I = J
J = J + J
ELSE
J = IR + 1
END IF
GOTO 20
END IF
RA(I) = RRA
GOTO 10
C
C 子程序SORT1结束
END
C
C
C
C=======================================================
SUBROUTINE INDEXX(N, ARRIN, INDX)
C
C
C 该子程序来自W.H. PRESS等人的“数字食谱”（"NUMERICAL RECIPES"）。
DIMENSION ARRIN(N), INDX(N)
C
DO 11 J = 1, N
INDX(J) = J
11 CONTINUE
L = (N / 2) + 1
IR = N
10 CONTINUE
IF (L .GT. 1) THEN
L = L - 1
INDXT = INDX(L)
Q = ARRIN(INDXT)
ELSE
INDXT = INDX(IR)
Q = ARRIN(INDXT)
INDX(IR) = INDX(1)
IR = IR - 1
IF (IR .EQ. 1) THEN
INDX(1) = INDXT
RETURN
END IF
END IF
I = L
J = L + L
20 IF (J .LE. IR) THEN
IF (J .LT. IR) THEN
IF (ARRIN(INDX(J)) .LT. ARRIN(INDX(J + 1))) J = J + 1
END IF
IF (Q .LT. ARRIN(INDX(J))) THEN
INDX(I) = INDX(J)
I = J
J = J + J
ELSE
J = IR + 1
END IF
GOTO 20
END IF
INDX(I) = INDXT
GOTO 10
C
C 子程序SORT1结束
END
C
C
C
C==============================================================
REAL FUNCTION RAN1(IDUM)
C
C
C 该子程序来自W.H. PRESS等人的“数字食谱”。
DIMENSION R(97)
PARAMETER (M1 = 259200, IA1 = 7141, IC1 = 54773, RM1 =
&3.8580247E-6)
PARAMETER (M2 = 134456, IA2 = 8121, IC2 = 28411, RM2 =
&7.4373773E-6)
PARAMETER (M3 = 243000, IA3 = 4561, IC3 = 51349)
SAVE
DATA IFF / 0 /
IF ((IDUM .LT. 0) .OR. (IFF .EQ. 0)) THEN
IFF = 1
IX1 = MOD(IC1 - IDUM,M1)
IX1 = MOD((IA1 * IX1) + IC1,M1)
IX2 = MOD(IX1,M2)
IX1 = MOD((IA1 * IX1) + IC1,M1)
IX3 = MOD(IX1,M3)
DO 11 J = 1, 97
IX1 = MOD((IA1 * IX1) + IC1,M1)
IX2 = MOD((IA2 * IX2) + IC2,M2)
R(J) = (DBLE(IX1) + (DBLE(IX2) * RM2)) * RM1
11 CONTINUE
IDUM = 1
END IF
IX1 = MOD((IA1 * IX1) + IC1,M1)
IX2 = MOD((IA2 * IX2) + IC2,M2)
IX3 = MOD((IA3 * IX3) + IC3,M3)
J = 1 + ((97 * IX3) / M3)
IF ((J .GT. 97) .OR. (J .LT. 1)) PAUSE
RAN1 = R(J)
R(J) = (DBLE(IX1) + (DBLE(IX2) * RM2)) * RM1
C
C 子程序RAN1结束
RETURN
END
C
C
C
C===============================================================
FUNCTION GASDEV(IDUM)
C
C
C 该子程序来自W.H. PRESS等人的“数字食谱”。
COMMON /GASBLK/ ISET
DATA ISET / 0 /
IF (ISET .EQ. 0) THEN
1 V1 = (2. * RAN1(IDUM)) - 1.
V2 = (2. * RAN1(IDUM)) - 1.
R = (V1 ** 2) + (V2 ** 2)
IF (R .GE. 1.) GOTO 1
FAC = SQRT(- ((2. * LOG(R)) / R))
GSET = V1 * FAC
GASDEV = V2 * FAC
ISET = 1
ELSE
GASDEV = GSET
ISET = 0
END IF
C
C 子程序GASDEV结束
RETURN
END
scein.for
SUBROUTINE SCEIN(NOPT,MAXN,KSTOP,PCENTO,ISEED,
& NGS,NPG,NPS,NSPL,MINGS,INIFLG,IPRINT)
C
C
C 此子程序读取并打印用于全局优化的
C 洗牌复形演化方法的输入变量
C ——版本2.1
C
C 作者：段青云 - 亚利桑那大学，1992年4月
C
C
CHARACTER*10 PCNTRL,DEFLT,USRSP
CHARACTER*4 REDUC,INITL,YSFLG,NOFLG
COMMON/IOPARS/ICNTRL,IOUT,IDAT,IWBAL,ISCE,IPE,IPC,IDET
C
DATA DEFLT/' DEFAULT '/
DATA USRSP/'USER SPEC.'/
DATA YSFLG/'YES '/
DATA NOFLG/'NO '/
C
IERROR = 0
IWARN = 0
WRITE(ISCE,700)
700 FORMAT(10X,'SHUFFLED COMPLEX EVOLUTION GLOBAL OPTIMIZATION',
& /,10X,46(1H=))
C
C
C 读入SCE控制参数
IDEFLT = 0
READ (ICNTRL,*)
READ (ICNTRL,800) MAXN,KSTOP,PCENTO,NGS,ISEED,IDEFLT
800 FORMAT(2I5,F5.2,3I5)
IF (ISEED .EQ. 0) ISEED = 1969
C
C
C 如果IDEFLT等于1，读入SCE控制参数
IF (IDEFLT .EQ. 1) THEN
READ(ICNTRL,810) NPG,NPS,NSPL,MINGS,INIFLG,IPRINT
810 FORMAT(6I5)
PCNTRL = USRSP
IF (NPG .EQ. 0) NPG = 2*NOPT + 1
IF (NPS .EQ. 0) NPS = NOPT + 1
IF (NSPL .EQ. 0) NSPL = NPG
IF (MINGS .EQ. 0) MINGS = NGS
ELSE
READ(ICNTRL,*)
PCNTRL = DEFLT
END IF
C
C
C 如果IDEFLT等于0，
C 则将SCE控制参数设置为默认值
IF (IDEFLT .EQ. 0) THEN
NPG = 2*NOPT + 1
NPS = NOPT + 1
NSPL = NPG
MINGS = NGS
INIFLG = 0
IPRINT = 0
END IF
C
C
C 检查SCE控制参数是否有效
IF (NGS .LT. 1 .OR. NGS .GE. 1320) THEN
WRITE(ISCE,900) NGS
900 FORMAT(//,1X,'**ERROR** NUMBER OF COMPLEXES IN INITIAL ',
* ' POPULATION ',I5,' IS NOT A VALID CHOICE')
IERROR = IERROR + 1
END IF
C
IF (KSTOP .LT. 0 .OR. KSTOP .GE. 10) THEN
WRITE(ISCE,901) KSTOP
901 FORMAT(//,1X,'**WARNING** THE NUMBER OF SHUFFLING LOOPS IN',
* ' WHICH THE CRITERION VALUE MUST CHANGE ',/,13X,'SHOULD BE',
* ' GREATER THAN 0 AND LESS THAN 10. ','KSTOP = ',I2,
* ' WAS SPECIFIED.'/,13X,'BUT KSTOP = 5 WILL BE USED INSTEAD.')
IWARN = IWARN + 1
KSTOP=5
END IF
C
IF (MINGS .LT. 1 .OR. MINGS .GT. NGS) THEN
WRITE(ISCE,902) MINGS
902 FORMAT(//,1X,'**WARNING** THE MINIMUM NUMBER OF COMPLEXES ',
* I2,' IS NOT A VALID CHOICE. SET IT TO DEFAULT')
IWARN = IWARN + 1
MINGS = NGS
END IF
C
IF (NPG .LT. 2 .OR. NPG .GT. 1320/MAX(NGS,1)) THEN
WRITE(ISCE,903) NPG
903 FORMAT(//,1X,'**WARNING** THE NUMBER OF POINTS IN A COMPLEX ',
* I4,' IS NOT A VALID CHOICE, SET IT TO DEFAULT')
IWARN = IWARN + 1
NPG = 2*NOPT+1
END IF
C
IF (NPS.LT.2 .OR. NPS.GT.NPG .OR. NPS.GT.50) THEN
WRITE(ISCE,904) NPS
904 FORMAT(//,1X,'**WARNING** THE NUMBER OF POINTS IN A SUB-',
* 'COMPLEX ',I4,' IS NOT A VALID CHOICE, SET IT TO DEFAULT')
IWARN = IWARN + 1
NPS = NOPT + 1
END IF
C
IF (NSPL .LT. 1) THEN
WRITE(ISCE,905) NSPL
905 FORMAT(//,1X,'**WARNING** THE NUMBER OF EVOLUTION STEPS ',
* 'TAKEN IN EACH COMPLEX BEFORE SHUFFLING ',I4,/,13X,
* 'IS NOT A VALID CHOICE, SET IT TO DEFAULT')
IWARN = IWARN + 1
NSPL = NPG
END IF
C
C 计算初始种群中的总点数
NPT = NGS * NPG
C
IF (NPT .GT. 1320) THEN
WRITE(ISCE,906) NPT
906 FORMAT(//,1X,'**WARNING** THE NUMBER OF POINTS IN INITIAL ',
* 'POPULATION ',I5,' EXCEED THE POPULATION LIMIT,',/,13X,
* 'SET NGS TO 2, AND NPG, NPS AND NSPL TO DEFAULTS')
IWARN = IWARN + 1
NGS = 2
NPG = 2*NOPT + 1
NPS = NOPT + 1
NSPL = NPG
END IF
C
C 打印错误和警告消息的总数
IF (IERROR .GE. 1) WRITE(ISCE,907) IERROR
907 FORMAT(//,1X,'*** TOTAL NUMBER OF ERROR MESSAGES IS ',I2)
C
IF (IWARN .GE. 1) WRITE(ISCE,908) IWARN
908 FORMAT(//,1X,'*** TOTAL NUMBER OF WARNING MESSAGES IS ',I2)
C
IF (MINGS .LT. NGS) THEN
REDUC = YSFLG
ELSE
REDUC = NOFLG
END IF
C
IF (INIFLG .NE. 0) THEN
INITL = YSFLG
ELSE
INITL = NOFLG
END IF
C
C
C 打印洗牌复形演化优化选项
104 WRITE(ISCE,910)
910 FORMAT(//,2X,'SCE CONTROL',5X,'MAX TRIALS',5X,
&'REQUIRED IMPROVEMENT',5X,'RANDOM',/,3X,'PARAMETER',8X,
&'ALLOWED',6X,'PERCENT',4X,'NO. LOOPS',6X,'SEED',/,
&2X,11(1H-),5X,10(1H-),5X,7(1H-),4X,9(1H-),5X,6(1H-))
C
PCENTA=PCENTO*100.
WRITE(ISCE,912) PCNTRL,MAXN,PCENTA,KSTOP,ISEED
912 FORMAT(3X,A10,7X,I5,10X,F3.1,9X,I2,9X,I5)
WRITE(ISCE,914) NGS,NPG,NPT,NPS,NSPL
914 FORMAT(//,18X,'SCE ALGORITHM CONTROL PARAMETERS',/,18X,32(1H=),
&//,2X,'NUMBER OF',5X,'POINTS PER',5X,'POINTS IN',6X,'POINTS PER',
&4X,'EVOL. STEPS',/,2X,'COMPLEXES',6X,'COMPLEX',6X,'INI. POPUL.',
&5X,'SUB-COMPLX',4X,'PER COMPLEX',/,2X,9(1H-),5X,10(1H-),4X,
&11(1H-),5X,10(1H-),4X,11(1H-),5X,/,2X,5(I5,10X))
WRITE(ISCE,915) REDUC,MINGS,INITL
915 FORMAT(//,15X,'COMPLX NO.',5X,'MIN COMPLEX',5X,'INI. POINT',/,
&15X,'REDUCTION',6X,'NO. ALLOWED',6X,'INCLUDED',/,
&15X,10(1H-),5X,11(1H-),5X,10(1H-),/,18X,A4,6X,I8,13X,A4)
IF (IERROR .GE. 1) THEN
WRITE(ISCE,922)
922 FORMAT(//,'*** THE OPTIMIZATION SEARCH IS NOT CONDUCTED BECAUSE',
& ' OF INPUT DATA ERROR ***')
STOP
END IF
C
C 子程序SCEIN结束
RETURN
END
参考代码
以下代码仅供参考，都无法成功运行。
Python版本
#-------------------------------------------------------------------------------
# Name: SCE_Python_shared version
# This is the implementation for the SCE algorithm,
# written by Q.Duan, 9/2004 - converted to python by Van Hoey S.2011
#-------------------------------------------------------------------------------
## Refer to paper:
## 'EFFECTIVE AND EFFICIENT GLOBAL OPTIMIZATION FOR CONCEPTUAL
## RAINFALL-RUNOFF MODELS' BY DUAN, Q., S. SOROOSHIAN, AND V.K. GUPTA,
## WATER RESOURCES RESEARCH, VOL 28(4), PP.1015-1031, 1992.
import random
import numpy as np
from ..test_functions import functn
import matplotlib.pyplot as plt
import pickle
import warnings
warnings.filterwarnings("ignore")
################################################################################
## Sampling called from SCE
################################################################################
def SampleInputMatrix(nrows,npars,bu,bl,iseed,distname='randomUniform'):
'''
Create input parameter matrix for nrows simulations,
for npars with bounds ub and lb (np.array from same size)
distname gives the initial sampling distribution (currently one for all parameters)
returns np.array
'''
np.random.seed(iseed)
x = np.zeros((nrows,npars))
bound = bu - bl
for i in range(nrows):
x[i,:] = bl + np.random.rand(1,npars)*bound
return x
def cceua(s,sf,bl,bu,icall,iseed,file=None):
# This is the subroutine for generating a new point in a simplex
nps,nopt=s.shape
n = nps
m = nopt
alpha = 1.0
beta = 0.5
# Assign the best and worst points:
sb=s[0,:]
fb=sf[0]
sw=s[-1,:]
fw=sf[-1]
if file != None:
model = open(file, 'rb')
functn = pickle.load(model)
model.close()
# Compute the centroid of the simplex excluding the worst point:
ce= np.mean(s[:-1,:],axis=0)
# Attempt a reflection point
snew = SampleInputMatrix(1,nopt,bu,bl,iseed,distname='randomUniform')
snew[0] = ce + alpha*(ce-sw)
# Check if is outside the bounds:
ibound=0
s1=snew[0]-bl
idx=(s1<0).nonzero()
if idx[0].size <> 0:
ibound=1
s1=bu-snew[0]
idx=(s1<0).nonzero()
if idx[0].size <> 0:
ibound=2
if ibound >= 1:
snew = SampleInputMatrix(1,nopt,bu,bl,iseed,distname='randomUniform')
global functn
fnew = functn.predict(snew)[0]
icall += 1
# Reflection failed; now attempt a contraction point:
if fnew > fw:
snew[0] = sw + beta * (ce-sw)
fnew = functn.predict(snew)[0]
icall += 1
# Both reflection and contraction have failed, attempt a random point;
if fnew > fw:
snew = SampleInputMatrix(1,nopt,bu,bl,iseed,distname='randomUniform')
fnew = functn.predict(snew)[0]
icall += 1
# END OF CCE
return snew[0],fnew,icall
def sceua(bl,bu,pf,ngs,file=None,plot=True):
# This is the subroutine implementing the SCE algorithm,
# written by Q.Duan, 9/2004 - converted to python by Van Hoey S.2011
# Initialize SCE parameters:
nopt=len(bl)
npg=2*nopt+1
nps=nopt+1
nspl=npg
npt=npg*ngs
bound=bu-bl
maxn=3000
kstop=10
pcento=0.1
peps=0.001
iseed=0
iniflg=0
ngs=2
if file != None:
model = open(file, 'rb')
functn = pickle.load(model)
model.close()
# Create an initial population to fill array x(npt,nopt):
x = SampleInputMatrix(npt,nopt,bu,bl,iseed,distname='randomUniform')
if iniflg==1:
x[0,:]=x0
nloop=0
icall=0
xf=np.zeros(npt)
global functn
xf=functn.predict(x)
for i in range(npt):
icall += 1
f0=xf[0]
# Sort the population in order of increasing function values;
idx = np.argsort(xf)
xf = np.sort(xf)
x=x[idx,:]
# Record the best and worst points;
bestx=x[0,:]
bestf=xf[0]
worstx=x[-1,:]
worstf=xf[-1]
BESTF=bestf
BESTX=bestx
ICALL=icall
# Compute the standard deviation for each parameter
xnstd=np.std(x,axis=0)
# Computes the normalized geometric range of the parameters
gnrng=np.exp(np.mean(np.log((np.max(x,axis=0)-np.min(x,axis=0))/bound)))
print 'The Initial Loop: 0'
print ' BESTF: %f ' %bestf
print ' BESTX: '
print bestx
print ' WORSTF: %f ' %worstf
print ' WORSTX: '
print worstx
print ' '
# Check for convergency;
if icall >= maxn:
print '*** OPTIMIZATION SEARCH TERMINATED BECAUSE THE LIMIT'
print 'ON THE MAXIMUM NUMBER OF TRIALS '
print maxn
print 'HAS BEEN EXCEEDED. SEARCH WAS STOPPED AT TRIAL NUMBER:'
print icall
print 'OF THE INITIAL LOOP!'
if gnrng < peps:
print 'THE POPULATION HAS CONVERGED TO A PRESPECIFIED SMALL PARAMETER SPACE'
# Begin evolution loops:
nloop = 0
criter=[]
criter_change=1e+5
while icall<maxn and gnrng>peps and criter_change>pcento:
nloop+=1
# Loop on complexes (sub-populations);
for igs in range(ngs):
# Partition the population into complexes (sub-populations);
cx=np.zeros((npg,nopt))
cf=np.zeros((npg))
k1=np.array(range(npg))
k2=k1*ngs+igs
cx[k1,:] = x[k2,:]
cf[k1] = xf[k2]
# Evolve sub-population igs for nspl steps:
for loop in range(nspl):
# Select simplex by sampling the complex according to a linear
# probability distribution
lcs=np.array([0]*nps)
lcs[0] = 1
for k3 in range(1,nps):
for i in range(1000):
## lpos = 1 + int(np.floor(npg+0.5-np.sqrt((npg+0.5)**2 - npg*(npg+1)*random.random())))
lpos = int(np.floor(npg+0.5-np.sqrt((npg+0.5)**2 - npg*(npg+1)*random.random())))
## idx=find(lcs(1:k3-1)==lpos)
idx=(lcs[0:k3]==lpos).nonzero() #check of element al eens gekozen
if idx[0].size == 0:
break
lcs[k3] = lpos
lcs.sort()
# Construct the simplex:
s = np.zeros((nps,nopt))
s=cx[lcs,:]
sf = cf[lcs]
snew,fnew,icall=cceua(s,sf,bl,bu,icall,iseed,file)
# Replace the worst point in Simplex with the new point:
s[-1,:] = snew
sf[-1] = fnew
# Replace the simplex into the complex;
cx[lcs,:] = s
cf[lcs] = sf
# Sort the complex;
idx = np.argsort(cf)
cf = np.sort(cf)
cx=cx[idx,:]
# End of Inner Loop for Competitive Evolution of Simplexes
#end of Evolve sub-population igs for nspl steps:
# Replace the complex back into the population;
x[k2,:] = cx[k1,:]
xf[k2] = cf[k1]
# End of Loop on Complex Evolution;
# Shuffled the complexes;
idx = np.argsort(xf)
xf = np.sort(xf)
x=x[idx,:]
PX=x
PF=xf
# Record the best and worst points;
bestx=x[0,:]
bestf=xf[0]
worstx=x[-1,:]
worstf=xf[-1]
BESTX = np.append(BESTX,bestx, axis=0)
BESTF = np.append(BESTF,bestf)
ICALL = np.append(ICALL,icall)
# Compute the standard deviation for each parameter
xnstd=np.std(x,axis=0)
# Computes the normalized geometric range of the parameters
gnrng=np.exp(np.mean(np.log((np.max(x,axis=0)-np.min(x,axis=0))/bound)))
print 'Evolution Loop: %d - Trial - %d' %(nloop,icall)
print ' BESTF: %f ' %bestf
print ' BESTX: '
print bestx
print ' WORSTF: %f ' %worstf
print ' WORSTX: '
print worstx
print ' '
# Check for convergency;
if icall >= maxn:
print '*** OPTIMIZATION SEARCH TERMINATED BECAUSE THE LIMIT'
print 'ON THE MAXIMUM NUMBER OF TRIALS '
print maxn
print 'HAS BEEN EXCEEDED.'
if gnrng < peps:
print 'THE POPULATION HAS CONVERGED TO A PRESPECIFIED SMALL PARAMETER SPACE'
criter=np.append(criter,bestf)
if nloop >= kstop:
criter_change= np.abs(criter[nloop-1]-criter[nloop-kstop])*100
criter_change= criter_change/np.mean(np.abs(criter[nloop-kstop:nloop]))
if criter_change < pcento:
print 'THE BEST POINT HAS IMPROVED IN LAST %d LOOPS BY LESS THAN THE THRESHOLD %f' %(kstop,pcento)
print 'CONVERGENCY HAS ACHIEVED BASED ON OBJECTIVE FUNCTION CRITERIA!!!'
# End of the Outer Loops
print 'SEARCH WAS STOPPED AT TRIAL NUMBER: %d' %icall
print 'NORMALIZED GEOMETRIC RANGE = %f' %gnrng
print 'THE BEST POINT HAS IMPROVED IN LAST %d LOOPS BY %f' %(kstop,criter_change)
#reshape BESTX
BESTX=BESTX.reshape(BESTX.size/nopt,nopt)
if plot:
fig = plt.figure()
ax1 = plt.subplot(111)
l1 = ax1.plot(BESTX)
ax1.legend((l1),(pf['names']), shadow=True)
# plt.title('Trace of the different parameters')
plt.xlabel('Evolution Loop')
plt.ylabel('Normalized Parameters\' value')
fig = plt.figure()
ax2 = plt.subplot(111)
ax2.plot(BESTF)
plt.title('Trace of model value')
plt.xlabel('Evolution Loop')
plt.ylabel('Model value')
plt.show()
# END of Subroutine sceua
return bestx,bestf,BESTX,BESTF,ICALL
MATLAB版本
以下代码仅供参考，都无法成功运行。
%
% Copyright (c) 2015, Mostapha Kalami Heris & Yarpiz (www.yarpiz.com)
% All rights reserved. Please read the "LICENSE" file for license terms.
%
% Project Code: YPEA110
% Project Title: Implementation of Shuffled Complex Evolution (SCE-UA)
% Publisher: Yarpiz (www.yarpiz.com)
%
% Developer: Mostapha Kalami Heris (Member of Yarpiz Team)
%
% Cite as:
% Mostapha Kalami Heris, Shuffled Complex Evolution in MATLAB (URL: https://yarpiz.com/80/ypea110-shuffled-complex-evolution), Yarpiz, 2015.
%
% Contact Info: sm.kalami@gmail.com, info@yarpiz.com
%
clc;
clear;
close all;
%% Problem Definition
% Objective Function
CostFunction = @(x) Sphere(x);
nVar = 10; % Number of Unknown Variables
VarSize = [1 nVar]; % Unknown Variables Matrix Size
VarMin = -10; % Lower Bound of Unknown Variables
VarMax = 10; % Upper Bound of Unknown Variables
%% SCE-UA Parameters
MaxIt = 500; % Maximum Number of Iterations
nPopComplex = 10; % Complex Size
nPopComplex = max(nPopComplex, nVar+1); % Nelder-Mead Standard
nComplex = 5; % Number of Complexes
nPop = nComplex*nPopComplex; % Population Size
I = reshape(1:nPop, nComplex, []);
% CCE Parameters
cce_params.q = max(round(0.5*nPopComplex), 2); % Number of Parents
cce_params.alpha = 3; % Number of Offsprings
cce_params.beta = 5; % Maximum Number of Iterations
cce_params.CostFunction = CostFunction;
cce_params.VarMin = VarMin;
cce_params.VarMax = VarMax;
%% Initialization
% Empty Individual Template
empty_individual.Position = [];
empty_individual.Cost = [];
% Initialize Population Array
pop = repmat(empty_individual, nPop, 1);
% Initialize Population Members
for i = 1:nPop
pop(i).Position = unifrnd(VarMin, VarMax, VarSize);
pop(i).Cost = CostFunction(pop(i).Position);
end
% Sort Population
pop = SortPopulation(pop);
% Update Best Solution Ever Found
BestSol = pop(1);
% Initialize Best Costs Record Array
BestCosts = nan(MaxIt, 1);
%% SCE-UA Main Loop
for it = 1:MaxIt
% Initialize Complexes Array
Complex = cell(nComplex, 1);
% Form Complexes and Run CCE
for j = 1:nComplex
% Complex Formation
Complex{j} = pop(I(j, :));
% Run CCE
Complex{j} = RunCCE(Complex{j}, cce_params);
% Insert Updated Complex into Population
pop(I(j, :)) = Complex{j};
end
% Sort Population
pop = SortPopulation(pop);
% Update Best Solution Ever Found
BestSol = pop(1);
% Store Best Cost Ever Found
BestCosts(it) = BestSol.Cost;
% Show Iteration Information
disp(['Iteration ' num2str(it) ': Best Cost = ' num2str(BestCosts(it))]);
end
%% Results
figure;
%plot(BestCosts, 'LineWidth', 2);
semilogy(BestCosts, 'LineWidth', 2);
xlabel('Iteration');
ylabel('Best Cost');
grid on;
C++版本
以下代码仅供参考，都无法成功运行。
SCEUA.h：
/******************************************************************************
File : SCEUA.h
Author : L. Shawn Matott - Converted from Origincal SCEUA Fortran code
Copyright : 2009, L. Shawn Matott
The SCE-UA method is a general purpose global optimization
program. It was originally developed by Dr. Qingyun Duan as part
of his doctoral dissertation work at the Department of Hydrology
and Water Resources, University of Arizona, Tucson, AZ 85721, USA.
The dissertation is entitled "A Global Optimization Strategy for
Efficient and Effective Calibration of Hydrologic Models". The
program has since been modified to make it easier for use on
problems of users' interests. The algorithm has been described
in detail in an article entitled "Effective and Efficient Global
Optimization for Conceptual Rainfall-Runoff Models", Water Resources
Research, Vol 28(4), pp.1015-1031, 1992; and in an article entitled
"A Shuffled Complex Evolution Approach for Effective and Efficient
Global Minimization" by Q. Duan, V.K. Gupta and S. Sorooshian,
Journal of Optimization Theory and its Applications, Vol 76(3),
pp 501-521, 1993. A paper entitled "Optimal Use of the SCE-UA Global
Optimization Method for Calibrating Watershed Models", by Q. Duan,
S. Sorooshian, & V.K. Gupta, Journal of Hydrology, Vol.158, 265-284,
1994, discussed how to use the SCE-UA Method in an efficient and
effective manner.
Input Summary for the SCEUA algorithm (adapted from original Fortran-based
description):
==========================================================================
variable type description
MAXN integer Maximum number of trials allowed before
optimization is terminated. The purpose of
MAXN is to stop an optimization search before
too much computer time is expended. MAXN
should be set large enough so that optimization
is generally completed before MAXN trials are
performed. Recommended value is 10,000 (increase or
decrease as necessary).
---> this parameter is called m_Budget within Ostrich
KSTOP integer Number of shuffling loops in which the
criterion must improve by the specified
percentage or else optimization will be
terminated. Recommended value: 5.
---> this parameter is called m_Kstop within Ostrich
PCENTO double Percentage by which the criterion value must
change in the specified number of shuffling
loops or else optimization is terminated
(Use decimal equivalent: Percentage/100).
Recommended value: 0.01.
---> this parameter is called m_Pcento within Ostrich
NGS integer Number of complexes used for optimization
search. Minimum value is 1.
Recommended value is between 2 and 20 depending
on the number of parameters to be optimized and
on the degree of difficulty of the problem.
---> this parameter is called m_NumComplexes within Ostrich
ISEED integer Random seed used in optimization search. Enter
any integer number. Default value (=1969) is
assumed if this field is left blank.
Recommended value: any large integer.
---> this parameter is called m_Seed within Ostrich
IDEFLT integer Flag for setting the control variables of the
SCE-UA algorithm. Enter false or leave the field
blank for default setting. Enter true for user
specified setting.
If option true is chosen, user must specify alg.
parameters.
---> this parameter is called m_bUserConfig within Ostrich
NPG integer Number of points in each complex. NPG should
be greater than or equal to 2. The default
value is equal to (2 * number of optimized
parameters + 1).
---> this parameter is called m_PtsPerComplex within Ostrich
NPS integer Number of points in each sub-complex. NPS
should be greater than or equal to 2 and less
than NPG. The default value is equal to
(number of optimized parameters + 1).
---> this parameter is called m_PtsPerSubComplex within Ostrich
NSPL integer Number of evolution steps taken by each complex
before next shuffling. Default value is equal
to NPG.
---> this parameter is called m_NumEvoSteps within Ostrich
MINGS integer Minimum number of complexes required for
optimization search, if the number of complexes
is allowed to reduce as the optimization search
proceeds. The default value is equal to NGS.
---> this parameter is called m_MinComplexes within Ostrich
INIFLG bool Flag on whether to include an initial point in
the starting population. Enter true if the initial
point is to be included. The default value is equal to false.
---> this parameter is called m_bUseInitPt within Ostrich
IPRINT integer Print-out control flag. Enter '0' to print out
the best estimate of the global optimum at the
end of each shuffling loop. Enter '1' to print
out every point in the entire sample population
at the end of each shuffling loop. The default
value is equal to 0. Enter 2 to ignore this variable
and use conventional Ostrich output.
---> this parameter is called m_OutputMode within Ostrich
PARAMS double * Initial estimates of the parameters to be optimized.
---> this parameter is called m_pParams within Ostrich
LOWER double * Lower bounds of the parameters to be optimized.
---> this parameter is called m_pLower within Ostrich
UPPER double * Upper bounds of the parameters to be optimized.
---> this parameter is called m_pUpper within Ostrich
Version History
10-31-09 lsm Created
******************************************************************************/
#ifndef SCEUA_H
#define SCEUA_H
// parent class
#include "AlgorithmABC.h"
// forward declarations
class ModelABC;
class StatsClass;
/******************************************************************************
class SCEUA
******************************************************************************/
class SCEUA : public AlgorithmABC
{
public:
SCEUA(ModelABC * pModel);
~SCEUA(void) { DBG_PRINT("SCEUA::DTOR"); Destroy(); }
void Destroy(void);
void Optimize(void);
void Calibrate(void);
void WriteMetrics(FILE * pFile);
void WarmStart(void);
int GetCurrentIteration(void) { return m_CurIter; }
private :
void InitFromFile(IroncladString pFileName);
void scemain(void);
void scein(double * a, double * bl, double * bu, int nopt, int *maxn,
int *kstop, double *pcento, int *iseed, int *ngs, int *npg,
int *nps, int *nspl, int *mings, int *iniflg, int *iprint);
void sceua(double * a, double * bl, double * bu, int nopt, int maxn,
int kstop, double pcento, int iseed, int ngs, int npg,
int nps, int nspl, int mings, int iniflg, int iprint);
void cce(int nopt, int nps, double ** s, double * sf, double * bl,
double * bu, double * xnstd, int * icall, double maxn,
int * iseed);
void getpnt(int nopt, int idist, int * iseed, double * x, double * bl,
double * bu, double * std, double * xi);
void sort(int n, int m, double ** rb, double * ra);
void parstt(int npt, int nopt, double ** x, double * xnstd,
double * bound, double * gnrng, int * ipcnvg);
void sort(int n, int * ra);
void comp(int n, int npt, int ngs1, int ngs2, int npg,
double ** a, double * af, double ** b, double * bf);
void chkcst(int nopt, double * snew, double * bl,
double * bu, int * ibound);
StatusStruct m_pStatus;
double m_Best;
int m_Budget; //MAXN
int m_Kstop; //KSTOP
double m_Pcento; //PCENTO
double m_Peps; //peps
double m_fSaved;
int m_NumComplexes; //NGS
int m_Seed; //ISEED
int m_UserConfig; //IDEFLT
int m_PtsPerComplex; //NPG
int m_PtsPerSubComplex; //NPS
int m_NumEvoSteps; //NSPL
int m_MinComplexes; //MINGS
int m_UseInitPt; //INIFLG
int m_OutputMode; //IPRINT
int m_np; //number of parameters
int m_CurIter;
double * m_pParams; //PARAMS
double * m_pLower; //LOWER
double * m_pUpper; //UPPER
bool m_bUseInitPt;
ModelABC * m_pModel;
StatsClass * m_pStats;
}; /* end class SCEUA */
extern "C" {
void SCEUA_Program(int argc, StringType argv[]);
}
#endif /* SCEUA_H */
SCEUA.cpp：
/******************************************************************************
File : SCEUA.cpp
Author : L. Shawn Matott - Converted from Origincal S/******************************************************************************
File : SCEUA.h
Author : L. Shawn Matott - Converted from Origincal SCEUA Fortran code
Copyright : 2009, L. Shawn Matott
The SCE-UA method is a general purpose global optimization
program. It was originally developed by Dr. Qingyun Duan as part
of his doctoral dissertation work at the Department of Hydrology
and Water Resources, University of Arizona, Tucson, AZ 85721, USA.
The dissertation is entitled "A Global Optimization Strategy for
Efficient and Effective Calibration of Hydrologic Models". The
program has since been modified to make it easier for use on
problems of users' interests. The algorithm has been described
in detail in an article entitled "Effective and Efficient Global
Optimization for Conceptual Rainfall-Runoff Models", Water Resources
Research, Vol 28(4), pp.1015-1031, 1992; and in an article entitled
"A Shuffled Complex Evolution Approach for Effective and Efficient
Global Minimization" by Q. Duan, V.K. Gupta and S. Sorooshian,
Journal of Optimization Theory and its Applications, Vol 76(3),
pp 501-521, 1993. A paper entitled "Optimal Use of the SCE-UA Global
Optimization Method for Calibrating Watershed Models", by Q. Duan,
S. Sorooshian, & V.K. Gupta, Journal of Hydrology, Vol.158, 265-284,
1994, discussed how to use the SCE-UA Method in an efficient and
effective manner.
Input Summary for the SCEUA algorithm (adapted from original Fortran-based
description):
==========================================================================
variable type description
MAXN integer Maximum number of trials allowed before
optimization is terminated. The purpose of
MAXN is to stop an optimization search before
too much computer time is expended. MAXN
should be set large enough so that optimization
is generally completed before MAXN trials are
performed. Recommended value is 10,000 (increase or
decrease as necessary).
---> this parameter is called m_Budget within Ostrich
KSTOP integer Number of shuffling loops in which the
criterion must improve by the specified
percentage or else optimization will be
terminated. Recommended value: 5.
---> this parameter is called m_Kstop within Ostrich
PCENTO double Percentage by which the criterion value must
change in the specified number of shuffling
loops or else optimization is terminated
(Use decimal equivalent: Percentage/100).
Recommended value: 0.01.
---> this parameter is called m_Pcento within Ostrich
NGS integer Number of complexes used for optimization
search. Minimum value is 1.
Recommended value is between 2 and 20 depending
on the number of parameters to be optimized and
on the degree of difficulty of the problem.
---> this parameter is called m_NumComplexes within Ostrich
ISEED integer Random seed used in optimization search. Enter
any integer number. Default value (=1969) is
assumed if this field is left blank.
Recommended value: any large integer.
---> this parameter is called m_Seed within Ostrich
IDEFLT integer Flag for setting the control variables of the
SCE-UA algorithm. Enter false or leave the field
blank for default setting. Enter true for user
specified setting.
If option true is chosen, user must specify alg.
parameters.
---> this parameter is called m_bUserConfig within Ostrich
NPG integer Number of points in each complex. NPG should
be greater than or equal to 2. The default
value is equal to (2 * number of optimized
parameters + 1).
---> this parameter is called m_PtsPerComplex within Ostrich
NPS integer Number of points in each sub-complex. NPS
should be greater than or equal to 2 and less
than NPG. The default value is equal to
(number of optimized parameters + 1).
---> this parameter is called m_PtsPerSubComplex within Ostrich
NSPL integer Number of evolution steps taken by each complex
before next shuffling. Default value is equal
to NPG.
---> this parameter is called m_NumEvoSteps within Ostrich
MINGS integer Minimum number of complexes required for
optimization search, if the number of complexes
is allowed to reduce as the optimization search
proceeds. The default value is equal to NGS.
---> this parameter is called m_MinComplexes within Ostrich
INIFLG bool Flag on whether to include an initial point in
the starting population. Enter true if the initial
point is to be included. The default value is equal to false.
---> this parameter is called m_bUseInitPt within Ostrich
IPRINT integer Print-out control flag. Enter '0' to print out
the best estimate of the global optimum at the
end of each shuffling loop. Enter '1' to print
out every point in the entire sample population
at the end of each shuffling loop. The default
value is equal to 0. Enter 2 to ignore this variable
and use conventional Ostrich output.
---> this parameter is called m_OutputMode within Ostrich
PARAMS double * Initial estimates of the parameters to be optimized.
---> this parameter is called m_pParams within Ostrich
LOWER double * Lower bounds of the parameters to be optimized.
---> this parameter is called m_pLower within Ostrich
UPPER double * Upper bounds of the parameters to be optimized.
---> this parameter is called m_pUpper within Ostrich
Version History
10-31-09 lsm Created
******************************************************************************/
#ifndef SCEUA_H
#define SCEUA_H
// parent class
#include "AlgorithmABC.h"
// forward declarations
class ModelABC;
class StatsClass;
/******************************************************************************
class SCEUA
******************************************************************************/
class SCEUA : public AlgorithmABC
{
public:
SCEUA(ModelABC * pModel);
~SCEUA(void) { DBG_PRINT("SCEUA::DTOR"); Destroy(); }
void Destroy(void);
void Optimize(void);
void Calibrate(void);
void WriteMetrics(FILE * pFile);
void WarmStart(void);
int GetCurrentIteration(void) { return m_CurIter; }
private :
void InitFromFile(IroncladString pFileName);
void scemain(void);
void scein(double * a, double * bl, double * bu, int nopt, int *maxn,
int *kstop, double *pcento, int *iseed, int *ngs, int *npg,
int *nps, int *nspl, int *mings, int *iniflg, int *iprint);
void sceua(double * a, double * bl, double * bu, int nopt, int maxn,
int kstop, double pcento, int iseed, int ngs, int npg,
int nps, int nspl, int mings, int iniflg, int iprint);
void cce(int nopt, int nps, double ** s, double * sf, double * bl,
double * bu, double * xnstd, int * icall, double maxn,
int * iseed);
void getpnt(int nopt, int idist, int * iseed, double * x, double * bl,
double * bu, double * std, double * xi);
void sort(int n, int m, double ** rb, double * ra);
void parstt(int npt, int nopt, double ** x, double * xnstd,
double * bound, double * gnrng, int * ipcnvg);
void sort(int n, int * ra);
void comp(int n, int npt, int ngs1, int ngs2, int npg,
double ** a, double * af, double ** b, double * bf);
void chkcst(int nopt, double * snew, double * bl,
double * bu, int * ibound);
StatusStruct m_pStatus;
double m_Best;
int m_Budget; //MAXN
int m_Kstop; //KSTOP
double m_Pcento; //PCENTO
double m_Peps; //peps
double m_fSaved;
int m_NumComplexes; //NGS
int m_Seed; //ISEED
int m_UserConfig; //IDEFLT
int m_PtsPerComplex; //NPG
int m_PtsPerSubComplex; //NPS
int m_NumEvoSteps; //NSPL
int m_MinComplexes; //MINGS
int m_UseInitPt; //INIFLG
int m_OutputMode; //IPRINT
int m_np; //number of parameters
int m_CurIter;
double * m_pParams; //PARAMS
double * m_pLower; //LOWER
double * m_pUpper; //UPPER
bool m_bUseInitPt;
ModelABC * m_pModel;
StatsClass * m_pStats;
}; /* end class SCEUA */
extern "C" {
void SCEUA_Program(int argc, StringType argv[]);
}
#endif /* SCEUA_H */CEUA Fortran code
Copyright : 2009, L. Shawn Matott
The SCE-UA method is a general purpose global optimization
program - Shuffled Complex Evolution. It was originally developed
by Dr. Qingyun Duan as part of his doctoral dissertation work at
University of Arizona, Tucson, AZ 85721, USA.
The dissertation is entitled "A Global Optimization Strategy for
Efficient and Effective Calibration of Hydrologic Models". The
program has since been modified to make it easier for use on
problems of users' interests.
The algorithm has been described in detail in an article entitled :
"Effective and Efficient Global Optimization for Conceptual Rainfall-Runoff Models",
Water Resources Research, Vol 28(4), pp.1015-1031, 1992
And in an article entitled:
"A Shuffled Complex Evolution Approach for Effective and Efficient Global Minimization",
by Q. Duan, V.K. Gupta and S. Sorooshian, Journal of Optimization Theory and its
Applications, Vol 76(3), pp 501-521, 1993.
A paper entitled "Optimal Use of the SCE-UA Global Optimization Method for Calibrating Watershed Models",
by Q. Duan, S. Sorooshian, & V.K. Gupta, Journal of Hydrology, Vol.158, 265-284, 1994,
discussed how to use the SCE-UA Method in an efficient and effective manner.
Input Summary for the SCEUA algorithm (adapted from original Fortran-based
description):
==========================================================================
variable type description
MAXN integer Maximum number of trials allowed before
optimization is terminated. The purpose of
MAXN is to stop an optimization search before
too much computer time is expended. MAXN
should be set large enough so that optimization
is generally completed before MAXN trials are
performed. Recommended value is 10,000 (increase or
decrease as necessary).
---> this parameter is called m_Budget within Ostrich
KSTOP integer Number of shuffling loops in which the
criterion must improve by the specified
percentage or else optimization will be
terminated. Recommended value: 5.
---> this parameter is called m_Kstop within Ostrich
PCENTO double Percentage by which the criterion value must
change in the specified number of shuffling
loops or else optimization is terminated
(Use decimal equivalent: Percentage/100).
Recommended value: 0.01.
---> this parameter is called m_Pcento within Ostrich
NGS integer Number of complexes used for optimization
search. Minimum value is 1.
Recommended value is between 2 and 20 depending
on the number of parameters to be optimized and
on the degree of difficulty of the problem.
---> this parameter is called m_NumComplexes within Ostrich
ISEED integer Random seed used in optimization search. Enter
any integer number. Default value (=1969) is
assumed if this field is left blank.
Recommended value: any large integer.
---> this parameter is called m_Seed within Ostrich
IDEFLT integer Flag for setting the control variables of the
SCE-UA algorithm. Enter false or leave the field
blank for default setting. Enter true for user
specified setting.
If option true is chosen, user must specify alg.
parameters.
---> this parameter is called m_bUserConfig within Ostrich
NPG integer Number of points in each complex. NPG should
be greater than or equal to 2. The default
value is equal to (2 * number of optimized
parameters + 1).
---> this parameter is called m_PtsPerComplex within Ostrich
NPS integer Number of points in each sub-complex. NPS
should be greater than or equal to 2 and less
than NPG. The default value is equal to
(number of optimized parameters + 1).
---> this parameter is called m_PtsPerSubComplex within Ostrich
NSPL integer Number of evolution steps taken by each complex
before next shuffling. Default value is equal
to NPG.
---> this parameter is called m_NumEvoSteps within Ostrich
MINGS integer Minimum number of complexes required for
optimization search, if the number of complexes
is allowed to reduce as the optimization search
proceeds. The default value is equal to NGS.
---> this parameter is called m_MinComplexes within Ostrich
INIFLG integer Flag on whether to include an initial point in
the starting population. Enter true if the initial
point is to be included. The default value is equal to false.
---> this parameter is called m_bUseInitPt within Ostrich
IPRINT integer Print-out control flag. Enter '0' to print out
the best estimate of the global optimum at the
end of each shuffling loop. Enter '1' to print
out every point in the entire sample population
at the end of each shuffling loop. The default
value is equal to 0. Enter 2 to ignore this variable
and use conventional Ostrich output.
---> this parameter is called m_OutputMode within Ostrich
PARAMS double * Initial estimates of the parameters to be optimized.
---> this parameter is called m_pParams within Ostrich
LOWER double * Lower bounds of the parameters to be optimized.
---> this parameter is called m_pLower within Ostrich
UPPER double * Upper bounds of the parameters to be optimized.
---> this parameter is called m_pUpper within Ostrich
Version History
10-31-09 lsm Created
******************************************************************************/
#include <math.h>
#include <string.h>
#include "SCEUA.h"
#include "Model.h"
#include "ParameterGroup.h"
#include "ParameterABC.h"
#include "StatsClass.h"
#include "Utility.h"
#include "WriteUtility.h"
#include "Exception.h"
/******************************************************************************
WarmStart()
Read the best solution from a previous run.
******************************************************************************/
void SCEUA::WarmStart(void)
{
int np = m_pModel->GetParamGroupPtr()->GetNumParams();
double * pbest = new double[np+1];
int newcount = SimpleWarmStart(np, pbest);
m_pModel->GetParamGroupPtr()->WriteParams(pbest);
((Model *)m_pModel)->SetCounter(newcount);
delete [] pbest;
}/* end WarmStart() */
/******************************************************************************
CTOR
Registers the algorithm pointer and creates instances of member variables.
******************************************************************************/
SCEUA::SCEUA(ModelABC * pModel)
{
RegisterAlgPtr(this);
m_pModel = pModel;
m_pParams = NULL;
m_pUpper = NULL;
m_pLower = NULL;
m_bUseInitPt = false;
m_fSaved = NEARLY_HUGE;
IncCtorCount();
}/* end CTOR() */
/******************************************************************************
Destroy()
******************************************************************************/
void SCEUA::Destroy(void)
{
delete [] m_pParams;
delete [] m_pLower;
delete [] m_pUpper;
}/* end Destroy() */
/******************************************************************************
Calibrate()
Solve the Least-Squares minimization problem using SCEUA.
******************************************************************************/
void SCEUA::Calibrate(void)
{
int id;
char fileName[DEF_STR_SZ];
FILE * pFile;
NEW_PRINT("StatsClass", 1);
m_pStats = new StatsClass(m_pModel);
MEM_CHECK(m_pStats);
RegisterStatsPtr(m_pStats);
Optimize();
//compute statistics (variance and covariance)
m_pStats->CalcStats();
id = 0;
sprintf(fileName, "OstOutput%d.txt", id);
//write statistics of best parameter set to output file
pFile = fopen(fileName, "a");
m_pStats->WriteStats(pFile);
fclose(pFile);
//write statistics of best parameter set to output file
m_pStats->WriteStats(stdout);
} /* end Calibrate() */
/******************************************************************************
Optimize()
Minimize the objective function using SCEUA.
******************************************************************************/
void SCEUA::Optimize(void)
{
ParameterGroup * pGroup = m_pModel->GetParamGroupPtr();
InitFromFile(GetInFileName());
WriteSetup(m_pModel, "Shuffled Complex Evolution - University of Arizona");
m_CurIter = 0;
//write banner
WriteBanner(m_pModel, "gen best value ", "Pct. Complete");
scemain(); //main SCE implemenation, converted from FORTRAN
//place model at optimal prameter set
pGroup->WriteParams(m_pParams);
m_pModel->Execute();
WriteOptimal(m_pModel, m_Best);
m_pStatus.pct = 100.00;
m_pStatus.numRuns = m_pModel->GetCounter();
WriteStatus(&m_pStatus);
//write algorithm metrics
WriteAlgMetrics(this);
} /* end Optimize() */
/******************************************************************************
scemain()
THIS IS THE MAIN PROGRAM CALLING SUBROUTINES SCEIN AND SCEUA
******************************************************************************/
void SCEUA::scemain(void)
{
double * a, * bl, * bu;
int jseed[10] = {2,3,5,7,11,13,17,19,23,29};
a = m_pParams;
bl = m_pLower;
bu = m_pUpper;
if(m_OutputMode != 2)
printf(" ENTER THE MAIN PROGRAM --- \n");
int nopt, kstop, iseed, ngs, npg, nps, nspl, mings, iprint, maxn;
int iniflg;
double pcento;
nopt = m_np = m_pModel->GetParamGroupPtr()->GetNumParams();
scein(a,bl,bu,nopt,&maxn,&kstop,&pcento,&iseed,
&ngs,&npg,&nps,&nspl,&mings,&iniflg,&iprint);
int nrun = 1;
int i;
for(i = 0; i < nrun; i++)
{
if (nrun != 1) iseed = jseed[i];
if(m_OutputMode != 2)
printf("@ SCE-UA Run Number %d Random Seed Value %d\n",i,iseed);
sceua(a,bl,bu,nopt,maxn,kstop,pcento,iseed,
ngs,npg,nps,nspl,mings,iniflg,iprint);
}/* end for() */
}/* end scemain() */
/******************************************************************************
scein()
THIS SUBROUTINE READS AND PRINTS THE INPUT VARIABLES FOR
SHUFFLED COMPLEX EVOLUTION METHOD FOR GLOBAL OPTIMIZATION
-- Version 2.1
ADAPTED FROM FORTRAN CODE WRITTEN BY
QINGYUN DUAN - UNIVERSITY OF ARIZONA, APRIL 1992
******************************************************************************/
void SCEUA::scein
(
double * a,
double * bl,
double * bu,
int nopt,
int *maxn,
int *kstop,
double *pcento,
int *iseed,
int *ngs,
int *npg,
int *nps,
int *nspl,
int *mings,
int *iniflg,
int *iprint
)
{
char pcntrl[100],deflt[100],usrsp[100];
char reduc[40],initl[40],ysflg[40],noflg[40],**xname;
strcpy(deflt, "DEFAULT");
strcpy(usrsp, "USER SPEC.");
strcpy(ysflg, "YES");
strcpy(noflg, "NO");
xname = new char*[nopt];
for(int i = 0; i < nopt; i++)
{
xname[i] = new char[50];
strncpy(xname[i], m_pModel->GetParamGroupPtr()->GetParamPtr(i)->GetName(), 9);
}
if(m_OutputMode != 2) printf("ENTER THE SCEIN SUBROUTINE --- \n");
//INITIALIZE I/O VARIABLES
FILE * pIn = fopen("sce.in", "r");
FILE * pOut = fopen("sce.out", "w");
int ierror = 0;
int iwarn = 0;
if(m_OutputMode != 2)
{
fprintf(pOut, "\
SHUFFLED COMPLEX EVOLUTION GLOBAL OPTIMIZATION\n\
==============================================\n\n\n");
}
// READ THE SCE CONTROL PARAMETERS
int ideflt = 0;
char line[160];
fgets(line, 1600, pIn);
sscanf(line, "%d %d %lf %d %d %d", maxn, kstop, pcento, ngs, iseed, &ideflt);
if (*iseed == 0) *iseed = 1969;
//IF ideflt IS EQUAL TO 1, READ THE SCE CONTROL PARAMETERS
if (ideflt == 1)
{
fgets(line, 1600, pIn);
sscanf(line, "%d %d %d %d %d %d", npg, nps, nspl, mings, iniflg, iprint);
strcpy(pcntrl, usrsp);
}
else
{
strcpy(pcntrl, deflt);
}
//READ THE INITIAL PARAMETER VALUES AND THE PARAMETER BOUNDS
int iopt;
for(iopt = 0; iopt < nopt; iopt++)
{
fgets(line, 1600, pIn);
sscanf(line,"%lf %lf %lf", &(a[iopt]), &(bl[iopt]), &(bu[iopt]));
}
//IF ideflt IS EQUAL TO 0, SET THE SCE CONTROL PARAMETERS TO THE DEFAULT VALUES
if (ideflt == 0)
{
*npg = 2*nopt + 1;
*nps = nopt + 1;
*nspl = *npg;
*mings = *ngs;
*iniflg = 0;
*iprint = 0;
}/* end if() */
//CHECK IF THE SCE CONTROL PARAMETERS ARE VALID
if ((*ngs < 1) || (*ngs >= 1320))
{
fprintf(pOut,
"**ERROR** NUMBER OF COMPLEXES IN INITIAL POPULATION (%d) IS NOT A VALID CHOICE\n",
*ngs);
ierror = ierror + 1;
}
if ((*kstop < 0) || (*kstop >= 20))
{
fprintf(pOut,
"**WARNING** THE NUMBER OF SHUFFLING LOOPS IN \
WHICH THE CRITERION VALUE MUST CHANGE SHOULD BE \
GREATER THAN 0 AND LESS THAN 10. kstop = %d WAS SPECIFIED. \
BUT kstop = 5 WILL BE USED INSTEAD.\n",
*kstop);
iwarn = iwarn + 1;
*kstop = 5;
}/* end if() */
if((*mings < 1) || (*mings > *ngs))
{
fprintf(pOut,
"**WARNING** THE MINIMUM NUMBER OF COMPLEXES (%d) \
IS NOT A VALID CHOICE. SET IT TO DEFAULT \n",
*mings);
iwarn = iwarn + 1;
*mings = *ngs;
}/* end if() */
if ((*npg < 2) || (*npg > 1320/MyMax(*ngs,1)))
{
fprintf(pOut,
"**WARNING** THE NUMBER OF POINTS IN A COMPLEX (%d) \
IS NOT A VALID CHOICE, SET IT TO DEFAULT\n",
*npg);
iwarn = iwarn + 1;
*npg = 2*nopt+1;
}/* end if() */
if ((*nps < 2) || (*nps > *npg) || (*nps > 50))
{
fprintf(pOut,
"**WARNING** THE NUMBER OF POINTS IN A SUB-COMPLEX (%d) \
IS NOT A VALID CHOICE, SET IT TO DEFAULT\n",
*nps);
iwarn = iwarn + 1;
*nps = nopt + 1;
}/* end if() */
if (*nspl < 1)
{
fprintf(pOut,
"**WARNING** THE NUMBER OF EVOLUTION STEPS \
TAKEN IN EACH COMPLEX BEFORE SHUFFLING (%d) \
IS NOT A VALID CHOICE, SET IT TO DEFAULT\n",
*nspl);
iwarn = iwarn + 1;
*nspl = *npg;
}/* end if() */
// COMPUTE THE TOTAL NUMBER OF POINTS IN INITIAL POPULATION
int npt;
npt = (*ngs) * (*npg); //npt = ngs * npg
if (npt > 1320)
{
fprintf(pOut,
"**WARNING** THE NUMBER OF POINTS IN INITIAL \
POPULATION (%d) EXCEED THE POPULATION LIMIT \
SET NGS TO 2, AND NPG, NPS AND NSPL TO DEFAULTS\n",
npt);
iwarn = iwarn + 1;
*ngs = 2;
*npg = 2*nopt + 1;
*nps = nopt + 1;
*nspl = *npg;
} /* end if() */
// PRINT OUT THE TOTAL NUMBER OF ERROR AND WARNING MESSAGES
if (ierror >= 1)
fprintf(pOut, "*** TOTAL NUMBER OF ERROR MESSAGES IS %d\n",ierror);
if (iwarn >= 1)
fprintf(pOut, "*** TOTAL NUMBER OF WARNING MESSAGES IS %d\n",iwarn);
if (*mings < *ngs)
strcpy(reduc, ysflg);
else
strcpy(reduc, noflg);
if (iniflg != 0)
strcpy(initl, ysflg);
else
strcpy(initl, noflg);
//PRINT SHUFFLED COMPLEX EVOLUTION OPTIMIZATION OPTIONS
fprintf(pOut,"\
SCE CONTROL MAX TRIALS REQUIRED IMPROVEMENT RANDOM\n\
PARAMETER ALLOWED PERCENT NO. LOOPS SEED\n\
----------- ---------- ------- --------- ------\n");
double pcenta;
pcenta=(*pcento)*100.;
fprintf(pOut," %-11s %-10d %7.2lf %-9d %-6d\n\n\n",
pcntrl, *maxn, pcenta, *kstop, *iseed);
fprintf(pOut,"\
SCE ALGORITHM CONTROL PARAMETERS\n\
================================\n\n\
NUMBER OF POINTS PER POINTS IN POINTS PER EVOL. STEPS\n\
COMPLEXES COMPLEX INI. POPUL. SUB-COMPLX PER COMPLEX\n\
--------- ---------- ----------- ---------- -----------\n");
fprintf(pOut," %-9d %-10d %-11d %-10d %-11d\n\n\n",
*ngs, *npg, npt, *nps, *nspl);
fprintf(pOut,"\
COMPLX NO. MIN COMPLEX INI. POINT\n\
REDUCTION NO. ALLOWED INCLUDED\n\
---------- ----------- ----------\n");
fprintf(pOut," %-10s %-11d %-10s\n\n\n",
reduc, *mings, initl);
fprintf(pOut,"\
INITIAL PARAMETER VALUES AND PARAMETER BOUNDS\n\
=============================================\n\n\
PARAMETER INITIAL VALUE LOWER BOUND UPPER BOUND\n\
--------- ------------- ----------- -----------\n");
for(int i = 0; i < nopt; i++)
{
fprintf(pOut, " %-9s %13.3lf %11.3lf %11.3lf\n",
xname[i], a[i], bl[i], bu[i]);
}/* end for() */
fprintf(pOut,"\n\n");
for(int i = 0; i < nopt; i++)
{
delete [] xname[i];
}
delete [] xname;
if (ierror >= 1)
{
fprintf(pOut,
"*** THE OPTIMIZATION SEARCH IS NOT CONDUCTED BECAUSE OF INPUT DATA ERROR ***\n");
fclose(pIn);
fclose(pOut);
ExitProgram(1);
}/* end if() */
fclose(pIn);
fclose(pOut);
}/* end scein() */
/******************************************************************************
sceua()
SHUFFLED COMPLEX EVOLUTION METHOD FOR GLOBAL OPTIMIZATION
-- Version 2.1
by QINGYUN DUAN (adpated to C++ by L. Shawn Matott)
DEPARTMENT OF HYDROLOGY & WATER RESOURCES
UNIVERSITY OF ARIZONA, TUCSON, AZ 85721
(602) 621-9360, email: duan@hwr.arizona.edu
WRITTEN IN OCTOBER 1990.
REVISED IN AUGUST 1991
REVISED IN APRIL 1992
STATEMENT BY AUTHOR:
--------------------
This general purpose global optimization program is developed at
the Department of Hydrology & Water Resources of the University
of Arizona. Further information regarding the SCE-UA method can
be obtained from Dr. Q. Duan, Dr. S. Sorooshian or Dr. V.K. Gupta
at the address and phone number listed above. We request all
users of this program make proper reference to the paper entitled
'Effective and Efficient Global Optimization for Conceptual
Rainfall-runoff Models' by Duan, Q., S. Sorooshian, and V.K. Gupta,
Water Resources Research, Vol 28(4), pp.1015-1031, 1992.
LIST OF INPUT ARGUEMENT VARIABLES
a(.) = initial parameter set
bl(.) = lower bound on parameters
bu(.) = upper bound on parameters
nopt = number of parameters to be optimized
LIST OF SCE ALGORITHMIC CONTROL PARAMETERS:
ngs = number of complexes in the initial population
npg = number of points in each complex
npt = total number of points in initial population (npt=ngs*npg)
nps = number of points in a sub-complex
nspl = number of evolution steps allowed for each complex before
complex shuffling
mings = minimum number of complexes required, if the number of
complexes is allowed to reduce as the optimization proceeds
iseed = initial random seed
iniflg = flag on whether to include the initial point in population
= 0, not included
= 1, included
iprint = flag for controlling print-out after each shuffling loop
= 0, print information on the best point of the population
= 1, print information on every point of the population
CONVERGENCE CHECK PARAMETERS
maxn = max no. of trials allowed before optimization is terminated
kstop = number of shuffling loops in which the criterion value must
chang by the given percentage before optimization is terminated
pcento = percentage by which the criterion value must change in
given number of shuffling loops
ipcnvg = flag indicating whether parameter convergence is reached
(i.e., check if gnrng is less than 0.001)
= 0, parameter convergence not satisfied
= 1, parameter convergence satisfied
LIST OF LOCAL VARIABLES
x(.,.) = coordinates of points in the population
xf(.) = function values of x(.,.)
xx(.) = coordinates of a single point in x
cx(.,.) = coordinates of points in a complex
cf(.) = function values of cx(.,.)
s(.,.) = coordinates of points in the current simplex
sf(.) = function values of s(.,.)
bestx(.) = best point at current shuffling loop
bestf = function value of bestx(.)
worstx(.) = worst point at current shuffling loop
worstf = function value of worstx(.)
xnstd(.) = standard deviation of parameters in the population
gnrng = normalized geometric mean of parameter ranges
lcs(.) = indices locating position of s(.,.) in x(.,.)
bound(.) = bound on ith variable being optimized
ngs1 = number of complexes in current population
ngs2 = number of complexes in last population
iseed1 = current random seed
criter(.) = vector containing the best criterion values of the last
10 shuffling loops
******************************************************************************/
void SCEUA::sceua
(
double * a,
double * bl,
double * bu,
int nopt,
int maxn,
int kstop,
double pcento,
int iseed,
int ngs,
int npg,
int nps,
int nspl,
int mings,
int iniflg,
int iprint
)
{
FILE * pOut = fopen("sce.out", "a");
//LOCAL ARRAYS
double ** x, * xx, * bestx, * worstx, * xf;
double ** s, * sf, ** cx, * cf;
double * xnstd, * bound, criter[20], * unit;
int * lcs;
//allocate memory
lcs = new int[nps];
sf = new double[nps];
xf = new double[ngs*npg];
cf = new double[npg];
x = new double * [ngs*npg];
cx = new double * [npg];
int i;
for(i = 0; i < ngs*npg; i++)
{
x[i] = new double[nopt];
}
for(i = 0; i < npg; i++)
{
cx[i] = new double[nopt];
}
xx = new double[nopt];
bestx = new double[nopt];
worstx = new double[nopt];
xnstd = new double[nopt];
bound = new double[nopt];
unit = new double[nopt];
s = new double *[nps];
for(i = 0; i < nps; i++)
{
s[i] = new double[nopt];
}
char **xname;
xname = new char*[nopt];
for(i = 0; i < nopt; i++)
{
xname[i] = new char[50];
strncpy(xname[i], m_pModel->GetParamGroupPtr()->GetParamPtr(i)->GetName(), 9);
}
if(m_OutputMode != 2) printf("ENTER THE SCEUA SUBROUTINE --- \n");
// INITIALIZE VARIABLES
int nloop, loop, igs, nopt1, nopt2, itmp;
nloop = 0;
loop = 0;
igs = 0;
nopt1 = 8;
if (nopt < 8) nopt1 = nopt;
nopt2 = 12;
if (nopt < 12) nopt2 = nopt;
//INITIALIZE RANDOM SEED TO A NEGATIVE INTEGER
int iseed1;
iseed1 = -abs(iseed);
//COMPUTE THE TOTAL NUMBER OF POINTS IN INITIAL POPUALTION
int npt, ngs1, npt1;
npt = ngs * npg;
ngs1 = ngs;
npt1 = npt;
fprintf(pOut, "\
==================================================\n\
ENTER THE SHUFFLED COMPLEX EVOLUTION GLOBAL SEARCH\n\
==================================================\n\n\n");
if(m_OutputMode != 2) printf(" *** Evolution Loop Number %d\n", nloop);
//COMPUTE THE BOUND FOR PARAMETERS BEING OPTIMIZED
int j;
for(j = 1; j <= nopt; j++)
{
bound[j-1] = bu[j-1] - bl[j-1];
unit[j-1] = 1.0;
}
//COMPUTE THE FUNCTION VALUE OF THE INITIAL POINT
double fa;
//handle warm start
if(m_pModel->CheckWarmStart() == true)
{
WarmStart();
m_pModel->GetParamGroupPtr()->ReadParams(a);
}
// handle parameter extraction
if(m_pModel->GetParamGroupPtr()->CheckExtraction() == true)
{
m_pModel->GetParamGroupPtr()->ReadParams(a);
}
m_pModel->GetParamGroupPtr()->WriteParams(a);
fa = m_pModel->Execute();
if (fa < m_fSaved)
{
m_fSaved = fa;
m_pModel->SaveBest(0);
}
//write initial config.
int nleft = m_Budget - m_pModel->GetCounter();
double eb = (double)(m_pModel->GetCounter())/(double)m_Budget; //elapsed budget
m_pStatus.curIter = 0;
m_pStatus.maxIter = m_Budget;
m_pStatus.pct = (float)100.00*((float)1.00-(float)nleft/(float)m_Budget);
m_pStatus.numRuns = m_pModel->GetCounter();
m_pModel->GetParamGroupPtr()->WriteParams(m_pParams);
WriteRecord(m_pModel, 0, fa, m_pStatus.pct);
m_CurIter++;
WriteStatus(&m_pStatus);
//PRINT THE INITIAL POINT AND ITS CRITERION VALUE
fprintf(pOut, "\
*** PRINT THE INITIAL POINT AND ITS CRITERION VALUE ***\n\n\
CRITERION ");
for(i = 0; i < nopt; i++)
{
fprintf(pOut, "%-9s ", xname[i]);
}
fprintf(pOut, "\n %8.0lf ", fa);
for(i = 0; i < nopt; i++)
{
fprintf(pOut, "%5.3lf ", a[i]);
}
fprintf(pOut, "\n\n\n");
// GENERATE AN INITIAL SET OF npt1 POINTS IN THE PARAMETER SPACE
if (iniflg == 1)
{
for(j = 0; j < nopt; j++)
x[0][j] = a[j];
xf[0] = fa;
WriteInnerEval(WRITE_SCE, npt, '.');
WriteInnerEval(1, npt, '.');
}/* end if() */
// ELSE, GENERATE A POINT RANDOMLY AND SET IT EQUAL TO x(1,.)
else
{
getpnt(nopt,1,&iseed1,xx,bl,bu,unit,bl);
eb = (double)(m_pModel->GetCounter())/(double)m_Budget;
for(j = 0; j < nopt; j++)
{
xx[j] = TelescopicCorrection(bl[j], bu[j], bestx[j], eb, xx[j]);
x[0][j] = xx[j];
}
WriteInnerEval(WRITE_SCE, npt, '.');
WriteInnerEval(1, npt, '.');
m_pModel->GetParamGroupPtr()->WriteParams(xx);
m_pModel->PerformParameterCorrections();
for(j = 0; j < nopt; j++)
{
xx[j] = m_pModel->GetParamGroupPtr()->GetParamPtr(j)->GetEstVal();
x[0][j] = xx[j];
}
xf[0] = m_pModel->Execute();
if (xf[0] < m_fSaved)
{
m_fSaved = xf[0];
m_pModel->SaveBest(0);
}
}
//use initial point if it's better than the randomly generated one
if(fa < xf[0])
{
fprintf(pOut, "THE INITIAL POINT IS BETTER THAN THE RANDOM STARTING POINT. USING IT INSTEAD.");
for(j = 0; j < nopt; j++)
x[0][j] = a[j];
xf[0] = fa;
}
int icall;
icall = 1;
if (icall >= maxn) goto label_9000;
//GENERATE npt1-1 RANDOM POINTS DISTRIBUTED UNIFORMLY IN THE PARAMETER
//SPACE, AND COMPUTE THE CORRESPONDING FUNCTION VALUES
label_restart:
for(i = 1; i < npt1; i++)
{
getpnt(nopt,1,&iseed1,xx,bl,bu,unit,bl);
eb = (double)(m_pModel->GetCounter())/(double)m_Budget;
for(j = 0; j < nopt; j++)
{
xx[j] = TelescopicCorrection(bl[j], bu[j], bestx[j], eb, xx[j]);
x[i][j] = xx[j];
}
WriteInnerEval(i+1, npt, '.');
m_pModel->GetParamGroupPtr()->WriteParams(xx);
m_pModel->PerformParameterCorrections();
for(j = 0; j < nopt; j++)
{
xx[j] = m_pModel->GetParamGroupPtr()->GetParamPtr(j)->GetEstVal();
x[i][j] = xx[j];
}
xf[i] = m_pModel->Execute();
if (xf[i] < m_fSaved)
{
m_fSaved = xf[i];
m_pModel->SaveBest(0);
}
icall = icall + 1;
if (icall >= maxn)
{
break;
}
}/* end for() */
WriteInnerEval(WRITE_ENDED, npt, '.');
// ARRANGE THE POINTS IN ORDER OF INCREASING FUNCTION VALUE
sort(npt1,nopt,x,xf);
for(j = 0; j < nopt; j++)
{
bestx[j] = x[0][j];
worstx[j] = x[npt1-1][j];
}
double bestf, worstf;
bestf = xf[0];
worstf = xf[npt1-1];
// COMPUTE THE PARAMETER RANGE FOR THE INITIAL POPULATION
double gnrng;
int ipcnvg;
parstt(npt1,nopt,x,xnstd,bound,&gnrng,&ipcnvg);
// PRINT THE RESULTS FOR THE INITIAL POPULATION
fprintf(pOut,"\
**** PRINT THE RESULTS OF THE SCE SEARCH ***\n\n\
LOOP TRIALS COMPLXS BEST F WORST F PAR RNG ");
for(i = 0; i < nopt; i++)
{
fprintf(pOut,"%-9s ", xname[i]);
}
fprintf(pOut,"\n");
fprintf(pOut," %4d %6d %7d %6.2lf %9.3E %8.3lf ",
nloop,icall,ngs1,bestf,worstf,gnrng);
for(i = 0; i < nopt; i++)
{
fprintf(pOut,"%6.3lf ", bestx[i]);
}
fprintf(pOut,"\n");
if (iprint == 1)
{
fprintf(pOut, "POPULATION AT LOOP (%d)\n", nloop);
for(i = 0; i < npt1; i++)
{
fprintf(pOut, "%8.3lf ", xf[i]);
for(j = 0; j < nopt; j++)
{
fprintf(pOut, "%8.3lf ", x[i][j]);
}
fprintf(pOut, "\n");
}/* end for() */
}/* end if() */
if (icall >= maxn) goto label_9000;
if (ipcnvg == 1) goto label_9200;
// BEGIN THE MAIN LOOP ----------------
label_1000:
nleft = m_Budget - m_pModel->GetCounter();
m_pStatus.curIter = nloop+1;
if(IsQuit() == true){ goto label_9999;}
if(nleft <= 0){ m_pStatus.pct = 100.00; goto label_9000;}
nloop = nloop + 1;
if(m_OutputMode != 2) printf(" *** Evolution Loop Number %d\n",nloop);
//BEGIN LOOP ON COMPLEXES
for(igs = 1; igs <= ngs1; igs++)
{
// ASSIGN POINTS INTO COMPLEXES
int k1, k2;
for(k1 = 1; k1 <= npg; k1++)
{
k2 = (k1-1) * ngs1 + igs;
for(j = 1; j <= nopt; j++)
{
cx[k1-1][j-1] = x[k2-1][j-1];
} //end for()
cf[k1-1] = xf[k2-1];
} //end for()
// BEGIN INNER LOOP - RANDOM SELECTION OF SUB-COMPLEXES ---------------
int tmp = 0;
WriteInnerEval(WRITE_SCE, m_NumEvoSteps, '.');
for(loop = 0; loop < nspl; loop++)
{
// CHOOSE A SUB-COMPLEX (nps points) ACCORDING TO A LINEAR
// PROBABILITY DISTRIBUTION
if (nps == npg)
{
int k;
for(k = 0; k < nps; k++)
{
lcs[k] = k;
} // end do
goto label_85;
} // end if
double myrand;
myrand = UniformRandom();
itmp = (int)(npg + 0.5 - sqrt(pow((npg+0.5),2.00) - npg*(npg+1.00)*myrand));
if(itmp >= npg) itmp = npg-1;
lcs[0] = itmp;
int k, lpos;
for(k = 2; k <= nps; k++)
{
label_60:
myrand = UniformRandom();
lpos = (int)(npg + 0.5 - sqrt(pow((npg+0.5),2.00) - npg*(npg+1.00)*myrand));
if(lpos >= npg) lpos = npg-1;
for(k1 = 1; k1 <= k-1; k1++)
{
if (lpos == lcs[k1-1]) goto label_60;
} // end do
lcs[k-1] = lpos;
} // end do
// ARRANGE THE SUB-COMPLEX IN ORDER OF INCEASING FUNCTION VALUE
sort(nps,lcs);
// CREATE THE SUB-COMPLEX ARRAYS
label_85:
for(k = 1; k <= nps; k++)
{
for(j = 1; j <= nopt; j++)
{
s[k-1][j-1] = cx[lcs[k-1]][j-1];
} // end do
sf[k-1] = cf[lcs[k-1]];
} // end do
// USE THE SUB-COMPLEX TO GENERATE NEW POINT(S)
cce(nopt,nps,s,sf,bl,bu,xnstd,&tmp,maxn,&iseed1);
// IF THE SUB-COMPLEX IS ACCEPTED, REPLACE THE NEW SUB-COMPLEX
// INTO THE COMPLEX
for(k = 1; k <= nps; k++)
{
for(j = 1; j <= nopt; j++)
{
cx[lcs[k-1]][j-1] = s[k-1][j-1];
} // end do
cf[lcs[k-1]] = sf[k-1];
} //end do
// SORT THE POINTS
sort(npg,nopt,cx,cf);
//IF MAXIMUM NUMBER OF RUNS EXCEEDED, BREAK OUT OF THE LOOP
if (icall >= maxn) break;
// END OF INNER LOOP ------------
} /* end for() */
WriteInnerEval(WRITE_ENDED, m_NumEvoSteps, '.');
icall += tmp;
// REPLACE THE NEW COMPLEX INTO ORIGINAL ARRAY x(.,.)
for(k1 = 1; k1 <= npg; k1++)
{
k2 = (k1-1) * ngs1 + igs;
for(j = 1; j <= nopt; j++)
{
x[k2-1][j-1] = cx[k1-1][j-1];
} // end do
xf[k2-1] = cf[k1-1];
} // end do
if (icall >= maxn) break;
//END LOOP ON COMPLEXES
} /* end for() */
// RE-SORT THE POINTS
sort(npt1,nopt,x,xf);
// RECORD THE BEST AND WORST POINTS
for(j = 0; j < nopt; j++)
{
m_pParams[j] = bestx[j] = x[0][j];
worstx[j] = x[npt1-1][j];
} //end do
m_Best = bestf = xf[0];
worstf = xf[npt1-1];
// TEST THE POPULATION FOR PARAMETER CONVERGENCE
parstt(npt1,nopt,x,xnstd,bound,&gnrng,&ipcnvg);
// PRINT THE RESULTS FOR CURRENT POPULATION
m_pModel->GetParamGroupPtr()->WriteParams(m_pParams);
nleft = m_Budget - m_pModel->GetCounter();
m_pStatus.pct = (float)100.00*((float)1.00-(float)nleft/(float)m_Budget);
m_pStatus.numRuns = m_pModel->GetCounter();
WriteStatus(&m_pStatus);
WriteRecord(m_pModel, nloop, m_Best, m_pStatus.pct);
m_CurIter++;
if((nloop%5) == 0)
{
fprintf(pOut,"\
LOOP TRIALS COMPLXS BEST F WORST F PAR RNG ");
for(i = 0; i < nopt; i++)
{
fprintf(pOut,"%-9s ", xname[i]);
}
fprintf(pOut,"\n");
}
fprintf(pOut," %4d %6d %7d %6.2lf %9.3E %8.3lf ",
nloop,icall,ngs1,bestf,worstf,gnrng);
for(i = 0; i < nopt; i++)
{
fprintf(pOut,"%6.3lf ", bestx[i]);
}
fprintf(pOut,"\n");
if (iprint == 1)
{
fprintf(pOut, "POPULATION AT LOOP (%d)\n", nloop);
for(i = 0; i < npt1; i++)
{
fprintf(pOut, "%8.3lf ", xf[i]);
for(j = 0; j < nopt; j++)
{
fprintf(pOut, "%8.3lf ", x[i][j]);
}
fprintf(pOut, "\n");
}/* end for() */
}/* end if() */
// TEST IF MAXIMUM NUMBER OF FUNCTION EVALUATIONS EXCEEDED
if (icall >= maxn) goto label_9000;
// COMPUTE THE COUNT ON SUCCESSIVE LOOPS W/O FUNCTION IMPROVEMENT
criter[19] = bestf;
if (nloop < (kstop+1)) goto label_132;
double denomi, timeou;
denomi = fabs(criter[19-kstop] + criter[19]) / 2.0;
timeou = fabs(criter[19-kstop] - criter[19]) / denomi;
if (timeou < pcento) goto label_9100;
label_132:
int l;
for(l=0; l < 19; l++)
{
criter[l] = criter[l+1];
} //end do
//IF POPULATION IS CONVERGED INTO A SUFFICIENTLY SMALL SPACE
if (ipcnvg == 1) goto label_9200;
//NONE OF THE STOPPING CRITERIA IS SATISFIED, CONTINUE SEARCH
//CHECK FOR COMPLEX NUMBER REDUCTION
int ngs2;
if (ngs1 > mings)
{
ngs2 = ngs1;
ngs1 = ngs1 - 1;
npt1 = ngs1 * npg;
comp(nopt,npt1,ngs1,ngs2,npg,x,xf,cx,cf);
} //end if
// END OF MAIN LOOP -----------
goto label_1000;
// SEARCH TERMINATED
label_9000:
fprintf(pOut, "\
*** OPTIMIZATION SEARCH TERMINATED BECAUSE THE\n\
LIMIT ON THE MAXIMUM NUMBER OF TRIALS (%d)\n\
WAS EXCEEDED. SEARCH WAS STOPPED AT %d SUB-COMPLEX\n\
OF COMPLEX %d IN SHUFFLING LOOP %d ***\n\n",
maxn, loop, igs, nloop);
goto label_9999;
label_9100:
fprintf(pOut, "\
*** OPTIMIZATION TERMINATED BECAUSE THE CRITERION\n\
VALUE HAS NOT CHANGED %5.2lf PERCENT IN %d\n\
SHUFFLING LOOPS ***\n\n", pcento*100.0,kstop);
goto label_9999;
label_9200:
fprintf(pOut, "\
*** OPTIMIZATION RESTARTED BECAUSE THE POPULATION HAS\n\
CONVERGED INTO %5.2lf PERCENT OF THE FEASIBLE SPACE ***\n\n",
gnrng*100.0);
goto label_restart;
label_9999:
// PRINT THE FINAL PARAMETER ESTIMATE AND ITS FUNCTION VALUE
fprintf(pOut, "\
*** PRINT THE FINAL PARAMETER ESTIMATE AND ITS CRITERION VALUE ***\n\n\
CRITERION ");
for(i = 0; i < nopt; i++)
{
fprintf(pOut,"%-9s ", xname[i]);
}
fprintf(pOut,"\n%6.3lf ", bestf);
for(i = 0; i < nopt; i++)
{
fprintf(pOut,"%6.3lf ", bestx[i]);
}
fprintf(pOut,"\n");
fclose(pOut);
/* clean up memory */
for(i = 0; i < nopt; i++)
{
delete [] xname[i];
}
delete [] xname;
for(i = 0; i < ngs*npg; i++)
{
delete [] x[i];
}
for(i = 0; i < npg; i++)
{
delete [] cx[i];
}
delete [] x;
delete [] cx;
delete [] xx;
delete [] bestx;
delete [] worstx;
delete [] xnstd;
delete [] bound;
delete [] unit;
for(i = 0; i < nps; i++)
{
delete [] s[i];
}
delete [] s;
delete [] sf;
delete [] lcs;
delete [] cf;
delete [] xf;
} /* end sceua() */
/******************************************************************************
cce()
ALGORITHM TO GENERATE A NEW POINT(S) FROM A SUB-COMPLEX
******************************************************************************/
void SCEUA::cce
(
int nopt,
int nps,
double ** s,
double * sf,
double * bl,
double * bu,
double * xnstd,
int * icall,
double maxn,
int * iseed
)
{
//SUB-COMPLEX VARIABLES
const double c1=0.8;
const double c2=0.4;
/* ----------------------------------------------------
LIST OF LOCAL VARIABLES
sb(.) = the best point of the simplex
sw(.) = the worst point of the simplex
w2(.) = the second worst point of the simplex
fw = function value of the worst point
ce(.) = the centroid of the simplex excluding wo
snew(.) = new point generated from the simplex
iviol = flag indicating if constraints are violated
= 1 , yes
= 0 , no
----------------------------------------------------- */
double * sw, * sb, *ce, *snew;
sw = new double[nopt];
sb = new double[nopt];
ce = new double[nopt];
snew = new double[nopt];
//EQUIVALENCE OF VARIABLES FOR READABILTY OF CODE
int n = nps;
int m = nopt;
double alpha = 1.0;
double beta = 0.5;
/* ---------------------------------------------------
IDENTIFY THE WORST POINT wo OF THE SUB-COMPLEX s
COMPUTE THE CENTROID ce OF THE REMAINING POINTS
COMPUTE step, THE VECTOR BETWEEN wo AND ce
IDENTIFY THE WORST FUNCTION VALUE fw
--------------------------------------------------- */
int i, j;
for(j = 0; j < m; j++)
{
sb[j] = s[0][j];
sw[j] = s[n-1][j];
ce[j] = 0.0;
for(i = 0; i < n-1; i++)
{
ce[j] = ce[j] + s[i][j];
} //end do
ce[j] = ce[j]/(double)(n-1);
} //end do
double fw;
fw = sf[n-1];
//COMPUTE THE NEW POINT snew
//FIRST TRY A REFLECTION STEP
for(j = 0; j < m; j++)
{
snew[j] = ce[j] + alpha * (ce[j] - sw[j]);
} //end do
//CHECK IF snew SATISFIES ALL CONSTRAINTS
int ibound;
chkcst(nopt,snew,bl,bu,&ibound);
/* ------------------------------------------------------------------
snew IS OUTSIDE THE BOUND,
CHOOSE A POINT AT RANDOM WITHIN FEASIBLE REGION ACCORDING TO
A NORMAL DISTRIBUTION WITH BEST POINT OF THE SUB-COMPLEX
AS MEAN AND STANDARD DEVIATION OF THE POPULATION AS STD
------------------------------------------------------------------ */
if (ibound >= 1) getpnt(nopt,2,iseed,snew,bl,bu,xnstd,sb);
//COMPUTE THE FUNCTION VALUE AT snew
double fnew;
WriteInnerEval(*icall+1, m_NumEvoSteps, '.');
double eb = (double)(m_pModel->GetCounter())/(double)m_Budget;
for(j = 0; j < m; j++) snew[j] = TelescopicCorrection(bl[j], bu[j], sb[j], eb, snew[j]);
m_pModel->GetParamGroupPtr()->WriteParams(snew);
m_pModel->PerformParameterCorrections();
for(j = 0; j < m; j++) snew[j] = m_pModel->GetParamGroupPtr()->GetParamPtr(j)->GetEstVal();
fnew = m_pModel->Execute();
if (fnew < m_fSaved)
{
m_fSaved = fnew;
m_pModel->SaveBest(0);
}
*icall = *icall + 1;
//COMPARE fnew WITH THE WORST FUNCTION VALUE fw
//fnew IS LESS THAN fw, ACCEPT THE NEW POINT snew AND RETURN
if (fnew <= fw) goto label_2000;
if (*icall >= maxn) goto label_3000;
//fnew IS GREATER THAN fw, SO TRY A CONTRACTION STEP
for(j = 0; j < m; j++)
{
snew[j] = ce[j] - beta * (ce[j] - sw[j]);
} //end do
//COMPUTE THE FUNCTION VALUE OF THE CONTRACTED POINT
eb = (double)(m_pModel->GetCounter())/(double)m_Budget;
for(j = 0; j < m; j++) snew[j] = TelescopicCorrection(bl[j], bu[j], sb[j], eb, snew[j]);
WriteInnerEval(*icall+1, m_NumEvoSteps, '.');
m_pModel->GetParamGroupPtr()->WriteParams(snew);
m_pModel->PerformParameterCorrections();
for(j = 0; j < m; j++) snew[j] = m_pModel->GetParamGroupPtr()->GetParamPtr(j)->GetEstVal();
fnew = m_pModel->Execute();
if (fnew < m_fSaved)
{
m_fSaved = fnew;
m_pModel->SaveBest(0);
}
*icall = *icall + 1;
//COMPARE fnew TO THE WORST VALUE fw
//IF fnew IS LESS THAN OR EQUAL TO fw, THEN ACCEPT THE POINT AND RETURN
if (fnew <= fw) goto label_2000;
if (*icall >= maxn) goto label_3000;
/* ---------------------------------------------------------------------
IF BOTH REFLECTION AND CONTRACTION FAIL, CHOOSE ANOTHER POINT
ACCORDING TO A NORMAL DISTRIBUTION WITH BEST POINT OF THE SUB-COMPLEX
AS MEAN AND STANDARD DEVIATION OF THE POPULATION AS STD
--------------------------------------------------------------------- */
getpnt(nopt,2,iseed,snew,bl,bu,xnstd,sb);
//COMPUTE THE FUNCTION VALUE AT THE RANDOM POINT
eb = (double)(m_pModel->GetCounter())/(double)m_Budget;
for(j = 0; j < m; j++) snew[j] = TelescopicCorrection(bl[j], bu[j], sb[j], eb, snew[j]);
WriteInnerEval(*icall+1, m_NumEvoSteps, '.');
m_pModel->GetParamGroupPtr()->WriteParams(snew);
m_pModel->PerformParameterCorrections();
for(j = 0; j < m; j++) snew[j] = m_pModel->GetParamGroupPtr()->GetParamPtr(j)->GetEstVal();
fnew = m_pModel->Execute();
if (fnew < m_fSaved)
{
m_fSaved = fnew;
m_pModel->SaveBest(0);
}
*icall = *icall + 1;
//REPLACE THE WORST POINT BY THE NEW POINT
label_2000:
for(j = 0; j < m; j++)
{
s[n-1][j] = snew[j];
} //end do
sf[n-1] = fnew;
label_3000:
//free up memory
delete [] sw;
delete [] sb;
delete [] ce;
delete [] snew;
} /* end cce() */
/******************************************************************************
getpnt()
This subroutine generates a new point within feasible region
x(.) = new point
xi(.) = focal point
bl(.) = lower bound
bu(.) = upper bound
std(.) = standard deviation of probability distribution
idist = probability flag
= 1 - uniform distribution
= 2 - Gaussian distribution
******************************************************************************/
void SCEUA::getpnt
(
int nopt,
int idist,
int * iseed,
double * x,
double * bl,
double * bu,
double * std,
double * xi
)
{
int ibound;
int j;
double myrand;
int icount = m_pModel->GetCounter();
label_1:
for(j = 0; j < nopt; j++) //1 do j=1, nopt
{
label_2:
if (idist == 1)
{
myrand = UniformRandom();
}
else if (idist == 2)
{
myrand = GaussRandom();
}
else
{
printf("unknown distribution!\n");
}
x[j] = xi[j] + std[j] * myrand * (bu[j] - bl[j]);
FILE * pFile = fopen("getpnt.txt", "a+");
fprintf(pFile, "%d\tx[%d]:%E\txi[%d]:%e\tstd[%d]:%E\tmyrand : %E\tbu[%d]:%E\tbl[%d]:%E\n",
icount, j+1, x[j], j+1, xi[j], j+1, std[j], myrand, j+1, bu[j], j+1, bl[j]);
fclose(pFile);
//Check explicit constraints
chkcst(1,&x[j],&bl[j],&bu[j],&ibound);
if (ibound >= 1)
{
goto label_2;
}
} //end do
//Check implicit constraints
chkcst(nopt,x,bl,bu,&ibound);
if (ibound >= 1)
{
goto label_1;
}
}/* end getpnt() */
/******************************************************************************
parstt()
SUBROUTINE CHECKING FOR PARAMETER CONVERGENCE
******************************************************************************/
void SCEUA::parstt
(
int npt,
int nopt,
double ** x,
double * xnstd,
double * bound,
double * gnrng,
int * ipcnvg
)
{
double * xmax, * xmin, * xmean;
xmax = new double[nopt];
xmin = new double[nopt];
xmean = new double[nopt];
const double delta = 1.0e-20;
const double peps = m_Peps;
//COMPUTE MAXIMUM, MINIMUM AND STANDARD DEVIATION OF PARAMETER VALUES
double gsum, xsum1, xsum2;
int i,k;
gsum = 0.0;
for(k = 0; k < nopt; k++)
{
xmax[k] = -1.0E+20;
xmin[k] = 1.0E+20;
xsum1 = 0.0;
xsum2 = 0.0;
for(i = 0; i < npt; i++)
{
xmax[k] = MyMax(x[i][k], xmax[k]);
xmin[k] = MyMin(x[i][k], xmin[k]);
xsum1 = xsum1 + x[i][k];
xsum2 = xsum2 + x[i][k]*x[i][k];
} //end do
xmean[k] = xsum1/(double)npt;
xnstd[k] = (xsum2 / (double)npt - xmean[k]*xmean[k]);
if (xnstd[k] <= delta) xnstd[k] = delta;
xnstd[k] = sqrt(xnstd[k]);
xnstd[k] = xnstd[k] / bound[k];
gsum = gsum + log(delta + (xmax[k]-xmin[k])/bound[k]);
} //end do
*gnrng = exp(gsum/(double)nopt);
//CHECK IF NORMALIZED STANDARD DEVIATION OF PARAMETER IS <= eps
*ipcnvg = 0;
if (*gnrng <= peps)
{
*ipcnvg = 1;
} //end if
delete [] xmax;
delete [] xmin;
delete [] xmean;
}/* end parstt() */
/******************************************************************************
comp()
THIS SUBROUTINE REDUCE INPUT MATRIX a(n,ngs2*npg) TO MATRIX
b(n,ngs1*npg) AND VECTOR af(ngs2*npg) TO VECTOR bf(ngs1*npg)
******************************************************************************/
void SCEUA::comp
(
int n,
int npt,
int ngs1,
int ngs2,
int npg,
double ** a,
double * af,
double ** b,
double * bf
)
{
int i, igs, ipg, k1, k2;
for(igs = 1; igs <= ngs1; igs++)
{
for(ipg = 1; ipg <= npg; ipg++)
{
k1=(ipg-1)*ngs2 + igs;
k2=(ipg-1)*ngs1 + igs;
while(k2 > npg) k2 -= npg;
for(i = 1; i <= n; i++)
{
b[k2-1][i-1] = a[k1-1][i-1];
} //end do
bf[k2-1] = af[k1-1];
} //end do
} //end do
int j, jj;
for(j = 0; j < npt; j++)
{
jj = j % npg;
for(i = 0; i < n; i++)
{
a[j][i] = b[jj][i];
} //end do
af[j] = bf[jj];
} //end do
} /* end comp() */
/******************************************************************************
sort()
Sort a complex of parameter sets in ascending order of cost function.
******************************************************************************/
void SCEUA::sort
(
int n,
int m,
double ** rb,
double * ra
)
{
//n = number of paramter sets
//m = number of parameters
//rb = list of parameter sets
//ra = list of costs
for(int p = 0; p < n; p++)
{
double fp = ra[p];
for(int i = p + 1; i < n; i++)
{
double fi = ra[i];
//swap if ith entry is better than pth
if(fi < fp)
{
ra[p] = fi;
ra[i] = fp;
for(int j = 0; j < m; j++)
{
double t = rb[i][j];
rb[i][j] = rb[p][j];
rb[p][j] = t;
}/* end for() */
}/* end if() */
}/* end for() */
}/* end for() */
} /* end sort() */
/******************************************************************************
sort()
Sort a list of integers in ascending order.
******************************************************************************/
void SCEUA::sort
(
int n,
int * ra
)
{
for(int p = 0; p < n; p++)
{
int fp = ra[p];
for(int i = p + 1; i < n; i++)
{
int fi = ra[i];
//swap if ith entry is better than pth
if(fi < fp)
{
ra[p] = fi;
ra[i] = fp;
}/* end if() */
}/* end for() */
}/* end for() */
} /* end sort() */
/******************************************************************************
chkcst()
This subroutine check if the trial point satisfies all
constraints.
ibound - violation indicator
= -1 initial value
= 0 no violation
= 1 violation
nopt = number of optimizing variables
ii = the ii'th variable of the arrays x, bl, and bu
******************************************************************************/
void SCEUA::chkcst
(
int nopt,
double * x,
double * bl,
double * bu,
int * ibound
)
{
*ibound = -1;
//Check if explicit constraints are violated
for(int ii=1; ii<=nopt; ii++)
{
if ((x[ii-1] < bl[ii-1]) || (x[ii-1] > bu[ii-1])) goto label10;
}
if (nopt == 1) goto label9;
// Check if implicit constraints are violated
// (no implicit constraints for this function)
//
// No constraints are violated
//
label9: *ibound = 0;
return;
// At least one of the constraints are violated
label10: *ibound = 1;
return;
}/* end chkcst() */
/******************************************************************************
InitFromFile()
Read configuration information from the given filename, then write the
configuration info. to the file "sce.in" (maintains compatibility with SCE
fortran implementation).
******************************************************************************/
void SCEUA::InitFromFile(IroncladString pFileName)
{
FILE * pFile;
int i;
char * line;
char tmp[DEF_STR_SZ], tmp2[DEF_STR_SZ];
//assign defaults
m_np = m_pModel->GetParamGroupPtr()->GetNumParams();
m_Budget = 10000; //MAXN
m_Kstop = 5; //KSTOP
m_Pcento = 0.01; //PCENTO
m_Peps = 1.0E-3; //peps
m_NumComplexes = (int)(sqrt((double)m_np)); //NGS
m_Seed = 1969; //ISEED
m_UserConfig = 1; //IDEFLT
m_PtsPerComplex = 2*m_np + 1; //NPG
m_PtsPerSubComplex = m_np+1; //NPS
m_NumEvoSteps = m_PtsPerComplex; //NSPL
m_MinComplexes = m_NumComplexes; //MINGS
m_UseInitPt = 1; //INIFLG
m_OutputMode = 2; //IPRINT
m_bUseInitPt = false; //INIFLG
//allocate initial parameter configuration
ParameterGroup * pGroup = m_pModel->GetParamGroupPtr();
m_np = pGroup->GetNumParams();
NEW_PRINT("double", m_np);
m_pParams = new double[m_np];
MEM_CHECK(m_pParams);
NEW_PRINT("double", m_np);
m_pLower = new double[m_np];
MEM_CHECK(m_pParams);
NEW_PRINT("double", m_np);
m_pUpper = new double[m_np];
MEM_CHECK(m_pParams);
for(i = 0; i < m_np; i++)
{
m_pParams[i] = pGroup->GetParamPtr(i)->GetEstVal();
m_pLower[i] = pGroup->GetParamPtr(i)->GetLwrBnd();
m_pUpper[i] = pGroup->GetParamPtr(i)->GetUprBnd();
}/* end for() */
//read in SCEUA configuration
pFile = fopen(pFileName, "r");
if(pFile == NULL)
{
//couldn't open file, use defaults and log the error.
LogError(ERR_FILE_IO, "Couldn't open SCEUA config. file. Using Defaults");
return;
}/* end if() */
if(CheckToken(pFile, "RandomSeed", GetInFileName()) == true)
{
fclose(pFile);
m_Seed = GetRandomSeed();
pFile = fopen(pFileName, "r");
}
rewind(pFile);
//make sure correct tokens are present
if(CheckToken(pFile, "BeginSCEUA", pFileName) == true)
{
FindToken(pFile, "EndSCEUA", pFileName);
rewind(pFile);
FindToken(pFile, "BeginSCEUA", pFileName);
line = GetNxtDataLine(pFile, pFileName);
while(strstr(line, "EndSCEUA") == NULL)
{
if(strstr(line, "Budget") != NULL)
{
sscanf(line, "%s %d", tmp, &m_Budget);
if(m_Budget < 100)
{
LogError(ERR_FILE_IO, "Invalid SCEUA budget. Defaulting to 100.");
m_Budget = 100;
}
}/*end if() */
else if(strstr(line, "LoopStagnationCriteria") != NULL)
{
sscanf(line, "%s %d", tmp, &m_Kstop);
}
else if(strstr(line, "PctChangeCriteria") != NULL)
{
sscanf(line, "%s %lf", tmp, &m_Pcento);
}
else if(strstr(line, "PopConvCriteria") != NULL)
{
sscanf(line, "%s %lf", tmp, &m_Peps);
}
else if(strstr(line, "NumComplexes") != NULL)
{
sscanf(line, "%s %d", tmp, &m_NumComplexes);
}
else if(strstr(line, "NumPointsPerComplex") != NULL)
{
sscanf(line, "%s %d", tmp, &m_PtsPerComplex);
}
else if(strstr(line, "NumPointsPerSubComplex") != NULL)
{
sscanf(line, "%s %d", tmp, &m_PtsPerSubComplex);
}
else if(strstr(line, "NumEvolutionSteps") != NULL)
{
sscanf(line, "%s %d", tmp, &m_NumEvoSteps);
}
else if(strstr(line, "MinNumOfComplexes") != NULL)
{
sscanf(line, "%s %d", tmp, &m_MinComplexes);
}
else if(strstr(line, "UseInitialPoint") != NULL)
{
sscanf(line, "%s %s", tmp, tmp2);
MyStrLwr(tmp2);
if(strcmp(tmp2, "yes") == 0) m_bUseInitPt = true;
}
else
{
sprintf(tmp, "Unknown token: %s", line);
LogError(ERR_FILE_IO, tmp);
}/* end else() */
line = GetNxtDataLine(pFile, pFileName);
} /* end while() */
}/* end if() */
fclose(pFile);
//create sce.in file
int iniflg = 0;
if(m_bUseInitPt) iniflg = 1;
FILE * pOut = fopen("sce.in", "w");
fprintf(pOut, "%d %d %lf %d %d 1\n",
m_Budget, m_Kstop, m_Pcento, m_NumComplexes, m_Seed);
fprintf(pOut, "%d %d %d %d %d 2\n",
m_PtsPerComplex, m_PtsPerSubComplex, m_NumEvoSteps, m_MinComplexes, iniflg);
for(i = 0; i < m_np; i++)
{
fprintf(pOut, "%E %E %E\n", m_pParams[i], m_pLower[i], m_pUpper[i]);
}
fclose(pOut);
} /* end InitFromFile() */
/******************************************************************************
WriteMetrics()
Write out algorithm metrics and setup.
******************************************************************************/
void SCEUA::WriteMetrics(FILE * pFile)
{
fprintf(pFile, "\nAlgorithm Metrics\n");
fprintf(pFile, "Algorithm : Shuffled Complex Evolution (SCE)\n");
fprintf(pFile, "Budget : %d\n", m_Budget);
fprintf(pFile, "Loop Stagnation Criteria : %d\n", m_Kstop);
fprintf(pFile, "Pct Change Criteria : %lf\n", m_Pcento);
fprintf(pFile, "Number of Complexes : %d\n", m_NumComplexes);
fprintf(pFile, "Points Per Complex : %d\n", m_PtsPerComplex);
fprintf(pFile, "Points Per Sub-Complex : %d\n", m_PtsPerSubComplex);
fprintf(pFile, "Num. of Evolution Steps : %d\n", m_NumEvoSteps);
fprintf(pFile, "Min. Num. of Complexes : %d\n", m_MinComplexes);
m_pModel->WriteMetrics(pFile);
}/* end WriteMetrics() */
/******************************************************************************
SCEUA_Program()
Calibrate or optimize the model using SCE.
******************************************************************************/
void SCEUA_Program(int argC, StringType argV[])
{
NEW_PRINT("Model", 1);
ModelABC * model = new Model();
NEW_PRINT("SCEUA", 1);
SCEUA * SCE = new SCEUA(model);
MEM_CHECK(SCE);
if(model->GetObjFuncId() == OBJ_FUNC_WSSE) { SCE->Calibrate(); }
else { SCE->Optimize(); }
delete SCE;
delete model;
} /* end SCEUA_Program() */
838
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
