【算法】02 SCE-UA简介及源代码
2024-06-17 19:30:31
基本原理
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++实现》下载
源码
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
CCHARACTER*4 XNAME(16)
C
C ARRAYS FROM THE INPUT DATADIMENSION A(16),BL(16),BU(16)
C
C LOCAL ARRAYSDIMENSION 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)
CCOMMON/IOPARS/ICNTRL,IOUT,IDAT,IWBAL,ISCE,IPE,IPC,IDET
CDATA XNAME /' X1',' X2',' X3',' X4',' X5',' X6',' X7',&' X8',' X9',' X10',' X11',' X12',' X13',' X14',' X15',' X16'/
C
C INITIALIZE VARIABLESWRITE(ISCE,400)400 FORMAT(//,2X,50(1H=),/,2X,'ENTER THE SHUFFLED COMPLEX EVOLUTION',& ' GLOBAL SEARCH',/,2X,50(1H=))WRITE(*,400)NLOOP = 0LOOP = 0IGS = 0
C
C INITIALIZE RANDOM SEED TO A NEGATIVE INTEGERISEED1 = -ABS(ISEED)
C
C COMPUTE THE TOTAL NUMBER OF POINTS IN INITIAL POPUALTIONNPT = NGS * NPGNGS1 = NGSNPT1 = NPT
C
C COMPUTE THE BOUND FOR PARAMETERS BEING OPTIMIZEDDO J = 1, NOPTBOUND(J) = BU(J) - BL(J)UNIT(J) = 1.0END DO
C
C COMPUTE THE FUNCTION VALUE OF THE INITIAL POINTFA = FUNCTN(NOPT,A)
C
C PRINT THE INITIAL POINT AND ITS CRITERION VALUEWRITE(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) THENDO J = 1, NOPTX(1,J) = A(J)END DOXF(1) = FA
C
C ELSE, GENERATE A POINT RANDOMLY AND SET IT EQUAL TO X(1,.)ELSECALL GETPNT(NOPT,1,ISEED1,XX,BL,BU,UNIT,BL)DO J=1, NOPTX(1,J) = XX(J)END DOXF(1) = FUNCTN(NOPT,XX)END IFICALL = 1IF (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 VALUESDO I = 2, NPT1CALL GETPNT(NOPT,1,ISEED1,XX,BL,BU,UNIT,BL)DO J = 1, NOPTX(I,J) = XX(J)END DOXF(I) = FUNCTN(NOPT,XX)ICALL = ICALL + 1IF (ICALL .GE. MAXN) THENNPT1 = IGO TO 45END IFEND DO
C
C ARRANGE THE POINTS IN ORDER OF INCREASING FUNCTION VALUE45 CALL SORT(NPT1,NOPT,X,XF)
C
C RECORD THE BEST AND WORST POINTSDO J = 1, NOPTBESTX(J) = X(1,J)WORSTX(J) = X(NPT1,J)END DOBESTF = XF(1)WORSTF = XF(NPT1)
C
C COMPUTE THE PARAMETER RANGE FOR THE INITIAL POPULATIONCALL PARSTT(NPT1,NOPT,X,XNSTD,BOUND,GNRNG,IPCNVG)
C
C PRINT THE RESULTS FOR THE INITIAL POPULATIONWRITE(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) THENWRITE(ISCE,650) NLOOPDO I = 1, NPT1WRITE(ISCE,660) XF(I),GNRNG,(X(I,J),J=1,NOPT)END DOEND IF
CIF (ICALL .GE. MAXN) GO TO 9000IF (IPCNVG .EQ. 1) GO TO 9200
C
C BEGIN THE MAIN LOOP ----------------1000 CONTINUENLOOP = NLOOP + 1
C
C BEGIN LOOP ON COMPLEXESDO 3000 IGS = 1, NGS1
C
C ASSIGN POINTS INTO COMPLEXESDO K1 = 1, NPGK2 = (K1-1) * NGS1 + IGSDO J = 1, NOPTCX(K1,J) = X(K2,J)END DOCF(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 DISTRIBUTIONIF (NPS .EQ. NPG) THENDO K = 1, NPSLCS(K) = KEND DOGO TO 85END IF
CRAND = RAN1(ISEED1)LCS(1) = 1 + INT(NPG + 0.5 - SQRT( (NPG+.5)**2 -& NPG * (NPG+1) * RAND ))DO K = 2, NPS60 RAND = RAN1(ISEED1)LPOS = 1 + INT(NPG + 0.5 - SQRT((NPG+.5)**2 -& NPG * (NPG+1) * RAND ))DO K1 = 1, K-1IF (LPOS .EQ. LCS(K1)) GO TO 60END DOLCS(K) = LPOSEND DO
C
C ARRANGE THE SUB-COMPLEX IN ORDER OF INCEASING FUNCTION VALUECALL SORT1(NPS,LCS)
C
C CREATE THE SUB-COMPLEX ARRAYS85 DO K = 1, NPSDO J = 1, NOPTS(K,J) = CX(LCS(K),J)END DOSF(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 COMPLEXDO K = 1, NPSDO J = 1, NOPTCX(LCS(K),J) = S(K,J)END DOCF(LCS(K)) = SF(K)END DO
C
C SORT THE POINTSCALL SORT(NPG,NOPT,CX,CF)
C
C IF MAXIMUM NUMBER OF RUNS EXCEEDED, BREAK OUT OF THE LOOPIF (ICALL .GE. MAXN) GO TO 2222
C
C END OF INNER LOOP ------------2000 CONTINUE2222 CONTINUE
C
C REPLACE THE NEW COMPLEX INTO ORIGINAL ARRAY X(.,.)DO K1 = 1, NPGK2 = (K1-1) * NGS1 + IGSDO J = 1, NOPTX(K2,J) = CX(K1,J)END DOXF(K2) = CF(K1)END DOIF (ICALL .GE. MAXN) GO TO 3333
C
C END LOOP ON COMPLEXES3000 CONTINUE
C
C RE-SORT THE POINTS3333 CALL SORT(NPT1,NOPT,X,XF)
C
C RECORD THE BEST AND WORST POINTSDO J = 1, NOPTBESTX(J) = X(1,J)WORSTX(J) = X(NPT1,J)END DOBESTF = XF(1)WORSTF = XF(NPT1)
C
C TEST THE POPULATION FOR PARAMETER CONVERGENCECALL PARSTT(NPT1,NOPT,X,XNSTD,BOUND,GNRNG,IPCNVG)
C
C PRINT THE RESULTS FOR CURRENT POPULATIONWRITE(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) THENWRITE(ISCE,650) NLOOPDO I = 1, NPT1WRITE(ISCE,660) XF(I),GNRNG,(X(I,J),J=1,NOPT)END DOEND IF
C
C TEST IF MAXIMUM NUMBER OF FUNCTION EVALUATIONS EXCEEDEDIF (ICALL .GE. MAXN) GO TO 9000
C
C COMPUTE THE COUNT ON SUCCESSIVE LOOPS W/O FUNCTION IMPROVEMENTCRITER(20) = BESTFIF (NLOOP .LT. (KSTOP+1)) GO TO 132DENOMI = ABS(CRITER(20-KSTOP) + CRITER(20)) / 2.TIMEOU = ABS(CRITER(20-KSTOP) - CRITER(20)) / DENOMIIF (TIMEOU .LT. PCENTO) GO TO 9100132 CONTINUEDO L = 1, 19CRITER(L) = CRITER(L+1)END DO
C
C IF POPULATION IS CONVERGED INTO A SUFFICIENTLY SMALL SPACEIF (IPCNVG .EQ. 1) GO TO 9200
C
C NONE OF THE STOPPING CRITERIA IS SATISFIED, CONTINUE SEARCH
C
C CHECK FOR COMPLEX NUMBER REDUCTIONIF (NGS1 .GT .MINGS) THENNGS2 = NGS1NGS1 = NGS1 - 1NPT1 = NGS1 * NPGCALL COMP(NOPT,NPT1,NGS1,NGS2,NPG,X,XF,CX,CF)END IF
C
C END OF MAIN LOOP -----------GO TO 1000
C
C SEARCH TERMINATED9000 CONTINUEWRITE(ISCE,800) MAXN,LOOP,IGS,NLOOPWRITE(*,800) MAXN,LOOP,IGS,NLOOPGO TO 99999100 CONTINUEWRITE(ISCE,810) PCENTO*100.,KSTOPWRITE(*,810) PCENTO*100.,KSTOPGO TO 99999200 WRITE(ISCE,820) GNRNG*100.WRITE(*,820) GNRNG*100.9999 CONTINUE
C
C PRINT THE FINAL PARAMETER ESTIMATE AND ITS FUNCTION VALUEWRITE(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 = BESTFDO J = 1, NOPTA(J) = BESTX(J)END DO
10000 CONTINUE
C
C END OF SUBROUTINE SCEUARETURN500 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 VARIABLESPARAMETER (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
CDIMENSION SW(16),SB(16),CE(16),SNEW(16)
C
C EQUIVALENCE OF VARIABLES FOR READABILTY OF CODEN = NPSM = NOPTALPHA = 1.0BETA = 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 FWDO J = 1, MSB(J) = S(1,J)SW(J) = S(N,J)CE(J) = 0.0DO I = 1, N-1CE(J) = CE(J) + S(I,J)END DOCE(J) = CE(J)/DBLE(N-1)END DOFW = SF(N)
C
C COMPUTE THE NEW POINT SNEW
C
C FIRST TRY A REFLECTION STEPDO J = 1, MSNEW(J) = CE(J) + ALPHA * (CE(J) - SW(J))END DO
C
C CHECK IF SNEW SATISFIES ALL CONSTRAINTSCALL 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 STDIF (IBOUND .GE. 1) CALL GETPNT(NOPT,2,ISEED,SNEW,BL,BU,XNSTD,SB)
C
C
C COMPUTE THE FUNCTION VALUE AT SNEWFNEW = 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 RETURNIF (FNEW .LE. FW) GO TO 2000IF (ICALL .GE. MAXN) GO TO 3000
C
C
C FNEW IS GREATER THAN FW, SO TRY A CONTRACTION STEPDO J = 1, MSNEW(J) = CE(J) - BETA * (CE(J) - SW(J))END DO
C
C COMPUTE THE FUNCTION VALUE OF THE CONTRACTED POINTFNEW = 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 RETURNIF (FNEW .LE. FW) GO TO 2000IF (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 STD1000 CALL GETPNT(NOPT,2,ISEED,SNEW,BL,BU,XNSTD,SB)
C
C COMPUTE THE FUNCTION VALUE AT THE RANDOM POINTFNEW = FUNCTN(NOPT,SNEW)ICALL = ICALL + 1
C
C
C REPLACE THE WORST POINT BY THE NEW POINT2000 CONTINUEDO J = 1, MS(N,J) = SNEW(J)END DOSF(N) = FNEW3000 CONTINUE
C
C END OF SUBROUTINE CCERETURNEND
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
CDIMENSION X(16),BL(16),BU(16),STD(16),XI(16)
C1 DO J=1, NOPT2 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 2END DO
C
C CHECK IMPLICIT CONSTRAINTS
C CALL CHKCST(NOPT,X,BL,BU,IBOUND)IF (IBOUND .GE. 1) GO TO 1
CRETURNEND
C
C
C
C===================================================================SUBROUTINE PARSTT(NPT,NOPT,X,XNSTD,BOUND,GNRNG,IPCNVG)
C
C SUBROUTINE CHECKING FOR PARAMETER CONVERGENCEDIMENSION 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 VALUESGSUM = 0.D0DO K = 1, NOPTXMAX(K) = -1.0D+20XMIN(K) = 1.0D+20XSUM1 = 0.D0XSUM2 = 0.D0DO I = 1, NPTXMAX(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 DOXMEAN(K) = XSUM1 / NPTXNSTD(K) = (XSUM2 / NPT - XMEAN(K)*XMEAN(K))IF (XNSTD(K) .LE. DELTA) XNSTD(K) = DELTAXNSTD(K) = SQRT(XNSTD(K))XNSTD(K) = XNSTD(K) / BOUND(K)GSUM = GSUM + LOG( DELTA + (XMAX(K)-XMIN(K))/BOUND(K) )END DOGNRNG = EXP(GSUM/NOPT)
C
C CHECK IF NORMALIZED STANDARD DEVIATION OF PARAMETER IS <= EPSIPCNVG = 0IF (GNRNG .LE. PEPS) THENIPCNVG = 1END IF
C
C END OF SUBROUTINE PARSTTRETURNEND
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, NGS1DO IPG=1, NPGK1=(IPG-1)*NGS2 + IGSK2=(IPG-1)*NGS1 + IGSDO I=1, NB(K2,I) = A(K1,I)END DOBF(K2) = AF(K1)END DOEND DO
CDO J=1, NPTDO I=1, NA(J,I) = B(J,I)END DOAF(J) = BF(J)END DO
C
C END OF SUBROUTINE COMPRETURNEND
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
CDIMENSION RA(2000),RB(2000,16),WK(2000,16),IWK(2000)
CCALL INDEXX(N, RA, IWK)DO 11 I = 1, NWK(I,1) = RA(I)11 CONTINUEDO 12 I = 1, NRA(I) = WK(IWK(I),1)12 CONTINUEDO 14 J = 1, MDO 13 I = 1, NWK(I,J) = RB(I,J)13 CONTINUE14 CONTINUEDO 16 J = 1, MDO 15 I = 1, NRB(I,J) = WK(IWK(I),J)15 CONTINUE16 CONTINUE
C
C END OF SUBROUTINE SORTRETURNEND
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
CDIMENSION RA(N)
CINTEGER RA, RRA
CL = (N / 2) + 1IR = N10 CONTINUEIF (L .GT. 1) THENL = L - 1RRA = RA(L)ELSERRA = RA(IR)RA(IR) = RA(1)IR = IR - 1IF (IR .EQ. 1) THENRA(1) = RRARETURNEND IFEND IFI = LJ = L + L20 IF (J .LE. IR) THENIF (J .LT. IR) THENIF (RA(J) .LT. RA(J + 1)) J = J + 1END IFIF (RRA .LT. RA(J)) THENRA(I) = RA(J)I = JJ = J + JELSEJ = IR + 1END IFGOTO 20END IFRA(I) = RRAGOTO 10
C
C END OF SUBROUTINE SORT1END
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)
CDO 11 J = 1, NINDX(J) = J11 CONTINUEL = (N / 2) + 1IR = N10 CONTINUEIF (L .GT. 1) THENL = L - 1INDXT = INDX(L)Q = ARRIN(INDXT)ELSEINDXT = INDX(IR)Q = ARRIN(INDXT)INDX(IR) = INDX(1)IR = IR - 1IF (IR .EQ. 1) THENINDX(1) = INDXTRETURNEND IFEND IFI = LJ = L + L20 IF (J .LE. IR) THENIF (J .LT. IR) THENIF (ARRIN(INDX(J)) .LT. ARRIN(INDX(J + 1))) J = J + 1END IFIF (Q .LT. ARRIN(INDX(J))) THENINDX(I) = INDX(J)I = JJ = J + JELSEJ = IR + 1END IFGOTO 20END IFINDX(I) = INDXTGOTO 10
C
C END OF SUBROUTINE INDEXXEND
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)SAVEDATA IFF / 0 /IF ((IDUM .LT. 0) .OR. (IFF .EQ. 0)) THENIFF = 1IX1 = 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, 97IX1 = MOD((IA1 * IX1) + IC1,M1)IX2 = MOD((IA2 * IX2) + IC2,M2)R(J) = (DBLE(IX1) + (DBLE(IX2) * RM2)) * RM111 CONTINUEIDUM = 1END IFIX1 = 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)) PAUSERAN1 = R(J)R(J) = (DBLE(IX1) + (DBLE(IX2) * RM2)) * RM1
C
C END OF SUBROUTINE RAN1RETURNEND
C
C
C
C===============================================================FUNCTION GASDEV(IDUM)
C
C
C THIS SUBROUTINE IS FROM "NUMERICAL RECIPES" BY PRESS ET AL.COMMON /GASBLK/ ISETDATA ISET / 0 /IF (ISET .EQ. 0) THEN1 V1 = (2. * RAN1(IDUM)) - 1.V2 = (2. * RAN1(IDUM)) - 1.R = (V1 ** 2) + (V2 ** 2)IF (R .GE. 1.) GOTO 1FAC = SQRT(- ((2. * LOG(R)) / R))GSET = V1 * FACGASDEV = V2 * FACISET = 1ELSEGASDEV = GSETISET = 0END IF
C
C END OF SUBROUTINE GASDEVRETURNEND
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
CCHARACTER*10 PCNTRL,DEFLT,USRSPCHARACTER*4 REDUC,INITL,YSFLG,NOFLGCOMMON/IOPARS/ICNTRL,IOUT,IDAT,IWBAL,ISCE,IPE,IPC,IDET
CDATA DEFLT/' DEFAULT '/DATA USRSP/'USER SPEC.'/DATA YSFLG/'YES '/DATA NOFLG/'NO '/
CIERROR = 0IWARN = 0WRITE(ISCE,700)700 FORMAT(10X,'SHUFFLED COMPLEX EVOLUTION GLOBAL OPTIMIZATION',& /,10X,46(1H=))
C
C
C READ THE SCE CONTROL PARAMETERSIDEFLT = 0READ (ICNTRL,*)READ (ICNTRL,800) MAXN,KSTOP,PCENTO,NGS,ISEED,IDEFLT800 FORMAT(2I5,F5.2,3I5)IF (ISEED .EQ. 0) ISEED = 1969
C
C
C IF IDEFLT IS EQUAL TO 1, READ THE SCE CONTROL PARAMETERSIF (IDEFLT .EQ. 1) THENREAD(ICNTRL,810) NPG,NPS,NSPL,MINGS,INIFLG,IPRINT810 FORMAT(6I5)PCNTRL = USRSPIF (NPG .EQ. 0) NPG = 2*NOPT + 1IF (NPS .EQ. 0) NPS = NOPT + 1IF (NSPL .EQ. 0) NSPL = NPGIF (MINGS .EQ. 0) MINGS = NGSELSEREAD(ICNTRL,*)PCNTRL = DEFLTEND IF
C
C
C IF IDEFLT IS EQUAL TO 0, SET THE SCE CONTROL PARAMETERS TO
C THE DEFAULT VALUESIF (IDEFLT .EQ. 0) THENNPG = 2*NOPT + 1NPS = NOPT + 1NSPL = NPGMINGS = NGSINIFLG = 0IPRINT = 0END IF
C
C
C CHECK IF THE SCE CONTROL PARAMETERS ARE VALIDIF (NGS .LT. 1 .OR. NGS .GE. 1320) THENWRITE(ISCE,900) NGS900 FORMAT(//,1X,'**ERROR** NUMBER OF COMPLEXES IN INITIAL ',* ' POPULATION ',I5,' IS NOT A VALID CHOICE')IERROR = IERROR + 1END IF
CIF (KSTOP .LT. 0 .OR. KSTOP .GE. 10) THENWRITE(ISCE,901) KSTOP901 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 + 1KSTOP=5END IF
CIF (MINGS .LT. 1 .OR. MINGS .GT. NGS) THENWRITE(ISCE,902) MINGS902 FORMAT(//,1X,'**WARNING** THE MINIMUM NUMBER OF COMPLEXES ',* I2,' IS NOT A VALID CHOICE. SET IT TO DEFAULT')IWARN = IWARN + 1MINGS = NGSEND IF
CIF (NPG .LT. 2 .OR. NPG .GT. 1320/MAX(NGS,1)) THENWRITE(ISCE,903) NPG903 FORMAT(//,1X,'**WARNING** THE NUMBER OF POINTS IN A COMPLEX ',* I4,' IS NOT A VALID CHOICE, SET IT TO DEFAULT')IWARN = IWARN + 1NPG = 2*NOPT+1END IF
CIF (NPS.LT.2 .OR. NPS.GT.NPG .OR. NPS.GT.50) THENWRITE(ISCE,904) NPS904 FORMAT(//,1X,'**WARNING** THE NUMBER OF POINTS IN A SUB-',* 'COMPLEX ',I4,' IS NOT A VALID CHOICE, SET IT TO DEFAULT')IWARN = IWARN + 1NPS = NOPT + 1END IF
CIF (NSPL .LT. 1) THENWRITE(ISCE,905) NSPL905 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 + 1NSPL = NPGEND IF
C
C COMPUTE THE TOTAL NUMBER OF POINTS IN INITIAL POPULATIONNPT = NGS * NPG
CIF (NPT .GT. 1320) THENWRITE(ISCE,906) NPT906 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 + 1NGS = 2NPG = 2*NOPT + 1NPS = NOPT + 1NSPL = NPGEND IF
C
C PRINT OUT THE TOTAL NUMBER OF ERROR AND WARNING MESSAGESIF (IERROR .GE. 1) WRITE(ISCE,907) IERROR907 FORMAT(//,1X,'*** TOTAL NUMBER OF ERROR MESSAGES IS ',I2)
CIF (IWARN .GE. 1) WRITE(ISCE,908) IWARN908 FORMAT(//,1X,'*** TOTAL NUMBER OF WARNING MESSAGES IS ',I2)
CIF (MINGS .LT. NGS) THENREDUC = YSFLGELSEREDUC = NOFLGEND IF
CIF (INIFLG .NE. 0) THENINITL = YSFLGELSEINITL = NOFLGEND IF
C
C
C PRINT SHUFFLED COMPLEX EVOLUTION OPTIMIZATION OPTIONS104 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-))
CPCENTA=PCENTO*100.WRITE(ISCE,912) PCNTRL,MAXN,PCENTA,KSTOP,ISEED912 FORMAT(3X,A10,7X,I5,10X,F3.1,9X,I2,9X,I5)WRITE(ISCE,914) NGS,NPG,NPT,NPS,NSPL914 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,INITL915 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) THENWRITE(ISCE,922)922 FORMAT(//,'*** THE OPTIMIZATION SEARCH IS NOT CONDUCTED BECAUSE',& ' OF INPUT DATA ERROR ***')STOPEND IF
C
C END OF SUBROUTINE SCEINRETURNEND翻译
以下代码仅供参考，都无法成功运行。
笔者根据源码翻译如下：
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 最佳标准值的向量
CCHARACTER*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)
CCOMMON/IOPARS/ICNTRL,IOUT,IDAT,IWBAL,ISCE,IPE,IPC,IDET
CDATA 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 = 0LOOP = 0IGS = 0
C
C 初始化随机种子为负整数ISEED1 = -ABS(ISEED)
C
C 计算在初始种群中的总点数NPT = NGS * NPGNGS1 = NGSNPT1 = NPT
C
C 计算待优化参数约束范围DO J = 1, NOPTBOUND(J) = BU(J) - BL(J)UNIT(J) = 1.0END 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) THENDO J = 1, NOPTX(1,J) = A(J)END DOXF(1) = FA
C
C 否则，随机生成一个点并设置其等于X(1,.)ELSECALL GETPNT(NOPT,1,ISEED1,XX,BL,BU,UNIT,BL)DO J=1, NOPTX(1,J) = XX(J)END DOXF(1) = FUNCTN(NOPT,XX)END IFICALL = 1IF (ICALL .GE. MAXN) GO TO 9000
C
C 生成在参数空间均匀分布的NPT1-1个随机点，
C 并计算对应的函数值DO I = 2, NPT1CALL GETPNT(NOPT,1,ISEED1,XX,BL,BU,UNIT,BL)DO J = 1, NOPTX(I,J) = XX(J)END DOXF(I) = FUNCTN(NOPT,XX)ICALL = ICALL + 1IF (ICALL .GE. MAXN) THENNPT1 = IGO TO 45END IFEND DO
C
C 按函数值递增的顺序排列点45 CALL SORT(NPT1,NOPT,X,XF)
C
C 记录最好和最差的点DO J = 1, NOPTBESTX(J) = X(1,J)WORSTX(J) = X(NPT1,J)END DOBESTF = 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) THENWRITE(ISCE,650) NLOOPDO I = 1, NPT1WRITE(ISCE,660) XF(I),GNRNG,(X(I,J),J=1,NOPT)END DOEND IF
CIF (ICALL .GE. MAXN) GO TO 9000IF (IPCNVG .EQ. 1) GO TO 9200
C
C 开始主循环 ----------------1000 CONTINUENLOOP = NLOOP + 1
C
C 在复形上开始循环DO 3000 IGS = 1, NGS1
C
C 为复形分配点DO K1 = 1, NPGK2 = (K1-1) * NGS1 + IGSDO J = 1, NOPTCX(K1,J) = X(K2,J)END DOCF(K1) = XF(K2)END DO
C
C 开始内循环 - 子复形的随机选择 ---------------DO 2000 LOOP = 1, NSPL
C
C 根据线性概率分布
C 选择子复形（NPS个点）IF (NPS .EQ. NPG) THENDO K = 1, NPSLCS(K) = KEND DOGO TO 85END IF
CRAND = RAN1(ISEED1)LCS(1) = 1 + INT(NPG + 0.5 - SQRT( (NPG+.5)**2 -& NPG * (NPG+1) * RAND ))DO K = 2, NPS60 RAND = RAN1(ISEED1)LPOS = 1 + INT(NPG + 0.5 - SQRT((NPG+.5)**2 -& NPG * (NPG+1) * RAND ))DO K1 = 1, K-1IF (LPOS .EQ. LCS(K1)) GO TO 60END DOLCS(K) = LPOSEND DO
C
C 按照函数值递增的顺序排列子复形CALL SORT1(NPS,LCS)
C
C 创建子复形数组85 DO K = 1, NPSDO J = 1, NOPTS(K,J) = CX(LCS(K),J)END DOSF(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, NPSDO J = 1, NOPTCX(LCS(K),J) = S(K,J)END DOCF(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 CONTINUE2222 CONTINUE
C
C 将新复形替换为原始数组 X(.,.)DO K1 = 1, NPGK2 = (K1-1) * NGS1 + IGSDO J = 1, NOPTX(K2,J) = CX(K1,J)END DOXF(K2) = CF(K1)END DOIF (ICALL .GE. MAXN) GO TO 3333
C
C 结束复形循环3000 CONTINUE
C
C 对点重新排序3333 CALL SORT(NPT1,NOPT,X,XF)
C
C 记录最好和最坏的点DO J = 1, NOPTBESTX(J) = X(1,J)WORSTX(J) = X(NPT1,J)END DOBESTF = 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) THENWRITE(ISCE,650) NLOOPDO I = 1, NPT1WRITE(ISCE,660) XF(I),GNRNG,(X(I,J),J=1,NOPT)END DOEND IF
C
C 测试是否超过了函数评估次数的最大值IF (ICALL .GE. MAXN) GO TO 9000
C
C 计算无函数改进的连续循环次数CRITER(20) = BESTFIF (NLOOP .LT. (KSTOP+1)) GO TO 132DENOMI = ABS(CRITER(20-KSTOP) + CRITER(20)) / 2.TIMEOU = ABS(CRITER(20-KSTOP) - CRITER(20)) / DENOMIIF (TIMEOU .LT. PCENTO) GO TO 9100132 CONTINUEDO L = 1, 19CRITER(L) = CRITER(L+1)END DO
C
C 如果种群聚集到一个足够小的空间IF (IPCNVG .EQ. 1) GO TO 9200
C
C 不满足任何停止标准，继续搜索
C
C 检查复形数量减少IF (NGS1 .GT .MINGS) THENNGS2 = NGS1NGS1 = NGS1 - 1NPT1 = NGS1 * NPGCALL COMP(NOPT,NPT1,NGS1,NGS2,NPG,X,XF,CX,CF)END IF
C
C 结束主循环 -----------GO TO 1000
C
C 搜索终止9000 CONTINUEWRITE(ISCE,800) MAXN,LOOP,IGS,NLOOPWRITE(*,800) MAXN,LOOP,IGS,NLOOPGO TO 99999100 CONTINUEWRITE(ISCE,810) PCENTO*100.,KSTOPWRITE(*,810) PCENTO*100.,KSTOPGO TO 99999200 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 = BESTFDO J = 1, NOPTA(J) = BESTX(J)END DO
10000 CONTINUE
C
C 子程序SCEUA结束RETURN500 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
CDIMENSION SW(16),SB(16),CE(16),SNEW(16)
C
C 代码可读性变量的等价性N = NPSM = NOPTALPHA = 1.0BETA = 0.5
C
C 确定子复形S的最坏点WO
C 计算剩余点的质心CE（CENTROID，简称CE）
C 计算步长，最坏点WO和质心CE之间的向量
C 确定最坏的函数值FW（The worst function value，简称FW）DO J = 1, MSB(J) = S(1,J)SW(J) = S(N,J)CE(J) = 0.0DO I = 1, N-1CE(J) = CE(J) + S(I,J)END DOCE(J) = CE(J)/DBLE(N-1)END DOFW = SF(N)
C
C 计算新点SNEW
C
C 首先尝试反射步骤DO J = 1, MSNEW(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 标准差作为STDIF (IBOUND .GE. 1) CALL GETPNT(NOPT,2,ISEED,SNEW,BL,BU,XNSTD,SB)
C
C
C 计算SNEW点的函数值FNEWFNEW = FUNCTN(NOPT,SNEW)ICALL = ICALL + 1
C
C 将FNEW与最差函数值FW进行比较
C
C FNEW小于FW，接受新点SNEW并返回IF (FNEW .LE. FW) GO TO 2000IF (ICALL .GE. MAXN) GO TO 3000
C
C
C FNEW大于FW，所以尝试收缩步骤DO J = 1, MSNEW(J) = CE(J) - BETA * (CE(J) - SW(J))END DO
C
C 计算收缩点的函数值FNEWFNEW = FUNCTN(NOPT,SNEW)ICALL = ICALL + 1
C
C 将FNEW与最差值FW进行比较
C 如果FNEW小于或等于FW，则接受该点并返回IF (FNEW .LE. FW) GO TO 2000IF (ICALL .GE. MAXN) GO TO 3000
C
C
C 如果反射和收缩都失败，
C 则根据正态分布选择另一个点，
C 子复形的最佳点作为种群的平均值并且标准偏差作为STD1000 CALL GETPNT(NOPT,2,ISEED,SNEW,BL,BU,XNSTD,SB)
C
C 计算随机点的函数值FNEW = FUNCTN(NOPT,SNEW)ICALL = ICALL + 1
C
C
C 用新点替换最坏点2000 CONTINUEDO J = 1, MS(N,J) = SNEW(J)END DOSF(N) = FNEW3000 CONTINUE
C
C 子程序CCE结束RETURNEND
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 - 高斯分布即正态分布
CDIMENSION X(16),BL(16),BU(16),STD(16),XI(16)
C1 DO J=1, NOPT2 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 2END DO
C
C 检查隐式约束
C CALL CHKCST(NOPT,X,BL,BU,IBOUND)IF (IBOUND .GE. 1) GO TO 1
CRETURNEND
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.D0DO K = 1, NOPTXMAX(K) = -1.0D+20XMIN(K) = 1.0D+20XSUM1 = 0.D0XSUM2 = 0.D0DO I = 1, NPTXMAX(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 DOXMEAN(K) = XSUM1 / NPTXNSTD(K) = (XSUM2 / NPT - XMEAN(K)*XMEAN(K))IF (XNSTD(K) .LE. DELTA) XNSTD(K) = DELTAXNSTD(K) = SQRT(XNSTD(K))XNSTD(K) = XNSTD(K) / BOUND(K)GSUM = GSUM + LOG( DELTA + (XMAX(K)-XMIN(K))/BOUND(K) )END DOGNRNG = EXP(GSUM/NOPT)
C
C 检查参数的标准化标准差是否 <= EPSIPCNVG = 0IF (GNRNG .LE. PEPS) THENIPCNVG = 1END IF
C
C 子程序PARSTT结束RETURNEND
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, NGS1DO IPG=1, NPGK1=(IPG-1)*NGS2 + IGSK2=(IPG-1)*NGS1 + IGSDO I=1, NB(K2,I) = A(K1,I)END DOBF(K2) = AF(K1)END DOEND DO
CDO J=1, NPTDO I=1, NA(J,I) = B(J,I)END DOAF(J) = BF(J)END DO
C
C 子程序COMP结束RETURNEND
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(.) = 局部变量
CDIMENSION RA(2000),RB(2000,16),WK(2000,16),IWK(2000)
CCALL INDEXX(N, RA, IWK)DO 11 I = 1, NWK(I,1) = RA(I)11 CONTINUEDO 12 I = 1, NRA(I) = WK(IWK(I),1)12 CONTINUEDO 14 J = 1, MDO 13 I = 1, NWK(I,J) = RB(I,J)13 CONTINUE14 CONTINUEDO 16 J = 1, MDO 15 I = 1, NRB(I,J) = WK(IWK(I),J)15 CONTINUE16 CONTINUE
C
C 子程序SORT结束RETURNEND
C
C
C===========================================================SUBROUTINE SORT1(N,RA)
C
C
C 改编自“数字食谱”（"NUMERICAL RECIPES"）的排序子程序
C 由 W.H. PRESS等人，PP. 231
C
C 变量列表
C RA(.) = 要排序的整数数组
CDIMENSION RA(N)
CINTEGER RA, RRA
CL = (N / 2) + 1IR = N10 CONTINUEIF (L .GT. 1) THENL = L - 1RRA = RA(L)ELSERRA = RA(IR)RA(IR) = RA(1)IR = IR - 1IF (IR .EQ. 1) THENRA(1) = RRARETURNEND IFEND IFI = LJ = L + L20 IF (J .LE. IR) THENIF (J .LT. IR) THENIF (RA(J) .LT. RA(J + 1)) J = J + 1END IFIF (RRA .LT. RA(J)) THENRA(I) = RA(J)I = JJ = J + JELSEJ = IR + 1END IFGOTO 20END IFRA(I) = RRAGOTO 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)
CDO 11 J = 1, NINDX(J) = J11 CONTINUEL = (N / 2) + 1IR = N10 CONTINUEIF (L .GT. 1) THENL = L - 1INDXT = INDX(L)Q = ARRIN(INDXT)ELSEINDXT = INDX(IR)Q = ARRIN(INDXT)INDX(IR) = INDX(1)IR = IR - 1IF (IR .EQ. 1) THENINDX(1) = INDXTRETURNEND IFEND IFI = LJ = L + L20 IF (J .LE. IR) THENIF (J .LT. IR) THENIF (ARRIN(INDX(J)) .LT. ARRIN(INDX(J + 1))) J = J + 1END IFIF (Q .LT. ARRIN(INDX(J))) THENINDX(I) = INDX(J)I = JJ = J + JELSEJ = IR + 1END IFGOTO 20END IFINDX(I) = INDXTGOTO 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)SAVEDATA IFF / 0 /IF ((IDUM .LT. 0) .OR. (IFF .EQ. 0)) THENIFF = 1IX1 = 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, 97IX1 = MOD((IA1 * IX1) + IC1,M1)IX2 = MOD((IA2 * IX2) + IC2,M2)R(J) = (DBLE(IX1) + (DBLE(IX2) * RM2)) * RM111 CONTINUEIDUM = 1END IFIX1 = 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)) PAUSERAN1 = R(J)R(J) = (DBLE(IX1) + (DBLE(IX2) * RM2)) * RM1
C
C 子程序RAN1结束RETURNEND
C
C
C
C===============================================================FUNCTION GASDEV(IDUM)
C
C
C 该子程序来自W.H. PRESS等人的“数字食谱”。COMMON /GASBLK/ ISETDATA ISET / 0 /IF (ISET .EQ. 0) THEN1 V1 = (2. * RAN1(IDUM)) - 1.V2 = (2. * RAN1(IDUM)) - 1.R = (V1 ** 2) + (V2 ** 2)IF (R .GE. 1.) GOTO 1FAC = SQRT(- ((2. * LOG(R)) / R))GSET = V1 * FACGASDEV = V2 * FACISET = 1ELSEGASDEV = GSETISET = 0END IF
C
C 子程序GASDEV结束RETURNENDscein.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
CCHARACTER*10 PCNTRL,DEFLT,USRSPCHARACTER*4 REDUC,INITL,YSFLG,NOFLGCOMMON/IOPARS/ICNTRL,IOUT,IDAT,IWBAL,ISCE,IPE,IPC,IDET
CDATA DEFLT/' DEFAULT '/DATA USRSP/'USER SPEC.'/DATA YSFLG/'YES '/DATA NOFLG/'NO '/
CIERROR = 0IWARN = 0WRITE(ISCE,700)700 FORMAT(10X,'SHUFFLED COMPLEX EVOLUTION GLOBAL OPTIMIZATION',& /,10X,46(1H=))
C
C
C 读入SCE控制参数IDEFLT = 0READ (ICNTRL,*)READ (ICNTRL,800) MAXN,KSTOP,PCENTO,NGS,ISEED,IDEFLT800 FORMAT(2I5,F5.2,3I5)IF (ISEED .EQ. 0) ISEED = 1969
C
C
C 如果IDEFLT等于1，读入SCE控制参数IF (IDEFLT .EQ. 1) THENREAD(ICNTRL,810) NPG,NPS,NSPL,MINGS,INIFLG,IPRINT810 FORMAT(6I5)PCNTRL = USRSPIF (NPG .EQ. 0) NPG = 2*NOPT + 1IF (NPS .EQ. 0) NPS = NOPT + 1IF (NSPL .EQ. 0) NSPL = NPGIF (MINGS .EQ. 0) MINGS = NGSELSEREAD(ICNTRL,*)PCNTRL = DEFLTEND IF
C
C
C 如果IDEFLT等于0，
C 则将SCE控制参数设置为默认值IF (IDEFLT .EQ. 0) THENNPG = 2*NOPT + 1NPS = NOPT + 1NSPL = NPGMINGS = NGSINIFLG = 0IPRINT = 0END IF
C
C
C 检查SCE控制参数是否有效IF (NGS .LT. 1 .OR. NGS .GE. 1320) THENWRITE(ISCE,900) NGS900 FORMAT(//,1X,'**ERROR** NUMBER OF COMPLEXES IN INITIAL ',* ' POPULATION ',I5,' IS NOT A VALID CHOICE')IERROR = IERROR + 1END IF
CIF (KSTOP .LT. 0 .OR. KSTOP .GE. 10) THENWRITE(ISCE,901) KSTOP901 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 + 1KSTOP=5END IF
CIF (MINGS .LT. 1 .OR. MINGS .GT. NGS) THENWRITE(ISCE,902) MINGS902 FORMAT(//,1X,'**WARNING** THE MINIMUM NUMBER OF COMPLEXES ',* I2,' IS NOT A VALID CHOICE. SET IT TO DEFAULT')IWARN = IWARN + 1MINGS = NGSEND IF
CIF (NPG .LT. 2 .OR. NPG .GT. 1320/MAX(NGS,1)) THENWRITE(ISCE,903) NPG903 FORMAT(//,1X,'**WARNING** THE NUMBER OF POINTS IN A COMPLEX ',* I4,' IS NOT A VALID CHOICE, SET IT TO DEFAULT')IWARN = IWARN + 1NPG = 2*NOPT+1END IF
CIF (NPS.LT.2 .OR. NPS.GT.NPG .OR. NPS.GT.50) THENWRITE(ISCE,904) NPS904 FORMAT(//,1X,'**WARNING** THE NUMBER OF POINTS IN A SUB-',* 'COMPLEX ',I4,' IS NOT A VALID CHOICE, SET IT TO DEFAULT')IWARN = IWARN + 1NPS = NOPT + 1END IF
CIF (NSPL .LT. 1) THENWRITE(ISCE,905) NSPL905 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 + 1NSPL = NPGEND IF
C
C 计算初始种群中的总点数NPT = NGS * NPG
CIF (NPT .GT. 1320) THENWRITE(ISCE,906) NPT906 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 + 1NGS = 2NPG = 2*NOPT + 1NPS = NOPT + 1NSPL = NPGEND IF
C
C 打印错误和警告消息的总数IF (IERROR .GE. 1) WRITE(ISCE,907) IERROR907 FORMAT(//,1X,'*** TOTAL NUMBER OF ERROR MESSAGES IS ',I2)
CIF (IWARN .GE. 1) WRITE(ISCE,908) IWARN908 FORMAT(//,1X,'*** TOTAL NUMBER OF WARNING MESSAGES IS ',I2)
CIF (MINGS .LT. NGS) THENREDUC = YSFLGELSEREDUC = NOFLGEND IF
CIF (INIFLG .NE. 0) THENINITL = YSFLGELSEINITL = NOFLGEND 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-))
CPCENTA=PCENTO*100.WRITE(ISCE,912) PCNTRL,MAXN,PCENTA,KSTOP,ISEED912 FORMAT(3X,A10,7X,I5,10X,F3.1,9X,I2,9X,I5)WRITE(ISCE,914) NGS,NPG,NPT,NPS,NSPL914 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,INITL915 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) THENWRITE(ISCE,922)922 FORMAT(//,'*** THE OPTIMIZATION SEARCH IS NOT CONDUCTED BECAUSE',& ' OF INPUT DATA ERROR ***')STOPEND IF
C
C 子程序SCEIN结束RETURNEND参考代码
以下代码仅供参考，都无法成功运行。
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 - blfor i in range(nrows):x[i,:] = bl + np.random.rand(1,npars)*boundreturn xdef cceua(s,sf,bl,bu,icall,iseed,file=None):# This is the subroutine for generating a new point in a simplexnps,nopt=s.shapen = npsm = noptalpha = 1.0beta = 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 pointsnew = SampleInputMatrix(1,nopt,bu,bl,iseed,distname='randomUniform')snew[0] = ce + alpha*(ce-sw)# Check if is outside the bounds:ibound=0s1=snew[0]-blidx=(s1<0).nonzero()if idx[0].size <> 0:ibound=1s1=bu-snew[0]idx=(s1<0).nonzero()if idx[0].size <> 0:ibound=2if ibound >= 1:snew = SampleInputMatrix(1,nopt,bu,bl,iseed,distname='randomUniform')global functnfnew = 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 CCEreturn snew[0],fnew,icalldef 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+1nps=nopt+1nspl=npgnpt=npg*ngsbound=bu-blmaxn=3000kstop=10pcento=0.1peps=0.001iseed=0iniflg=0ngs=2if 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,:]=x0nloop=0icall=0xf=np.zeros(npt)global functnxf=functn.predict(x)for i in range(npt):icall += 1f0=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=bestfBESTX=bestxICALL=icall# Compute the standard deviation for each parameterxnstd=np.std(x,axis=0)# Computes the normalized geometric range of the parametersgnrng=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 ' %bestfprint ' BESTX: 'print bestxprint ' WORSTF: %f ' %worstfprint ' WORSTX: 'print worstxprint ' '# Check for convergency;if icall >= maxn:print '*** OPTIMIZATION SEARCH TERMINATED BECAUSE THE LIMIT'print 'ON THE MAXIMUM NUMBER OF TRIALS 'print maxnprint 'HAS BEEN EXCEEDED. SEARCH WAS STOPPED AT TRIAL NUMBER:'print icallprint 'OF THE INITIAL LOOP!'if gnrng < peps:print 'THE POPULATION HAS CONVERGED TO A PRESPECIFIED SMALL PARAMETER SPACE'# Begin evolution loops:nloop = 0criter=[]criter_change=1e+5while 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+igscx[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 distributionlcs=np.array([0]*nps)lcs[0] = 1for 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 gekozenif idx[0].size == 0:breaklcs[k3] = lposlcs.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,:] = snewsf[-1] = fnew# Replace the simplex into the complex;cx[lcs,:] = scf[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=xPF=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 parameterxnstd=np.std(x,axis=0)# Computes the normalized geometric range of the parametersgnrng=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 ' %bestfprint ' BESTX: 'print bestxprint ' WORSTF: %f ' %worstfprint ' WORSTX: 'print worstxprint ' '# Check for convergency;if icall >= maxn:print '*** OPTIMIZATION SEARCH TERMINATED BECAUSE THE LIMIT'print 'ON THE MAXIMUM NUMBER OF TRIALS 'print maxnprint '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])*100criter_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 Loopsprint 'SEARCH WAS STOPPED AT TRIAL NUMBER: %d' %icallprint 'NORMALIZED GEOMETRIC RANGE = %f' %gnrngprint 'THE BEST POINT HAS IMPROVED IN LAST %d LOOPS BY %f' %(kstop,criter_change)#reshape BESTXBESTX=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 sceuareturn 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 SizeVarMin = -10; % Lower Bound of Unknown Variables
VarMax = 10; % Upper Bound of Unknown Variables%% SCE-UA ParametersMaxIt = 500; % Maximum Number of IterationsnPopComplex = 10; % Complex Size
nPopComplex = max(nPopComplex, nVar+1); % Nelder-Mead StandardnComplex = 5; % Number of Complexes
nPop = nComplex*nPopComplex; % Population SizeI = 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:nPoppop(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 Loopfor it = 1:MaxIt% Initialize Complexes ArrayComplex = cell(nComplex, 1);% Form Complexes and Run CCEfor j = 1:nComplex% Complex FormationComplex{j} = pop(I(j, :));% Run CCEComplex{j} = RunCCE(Complex{j}, cce_params);% Insert Updated Complex into Populationpop(I(j, :)) = Complex{j};end% Sort Populationpop = SortPopulation(pop);% Update Best Solution Ever FoundBestSol = pop(1);% Store Best Cost Ever FoundBestCosts(it) = BestSol.Cost;% Show Iteration Informationdisp(['Iteration ' num2str(it) ': Best Cost = ' num2str(BestCosts(it))]);end%% Resultsfigure;
%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 beforeoptimization is terminated. The purpose ofMAXN is to stop an optimization search beforetoo much computer time is expended. MAXNshould be set large enough so that optimizationis generally completed before MAXN trials areperformed. Recommended value is 10,000 (increase ordecrease as necessary).
---> this parameter is called m_Budget within Ostrich
KSTOP integer Number of shuffling loops in which the criterion must improve by the specifiedpercentage or else optimization will beterminated. Recommended value: 5.
---> this parameter is called m_Kstop within Ostrich
PCENTO double Percentage by which the criterion value mustchange 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 optimizationsearch. Minimum value is 1.Recommended value is between 2 and 20 dependingon the number of parameters to be optimized andon the degree of difficulty of the problem.
---> this parameter is called m_NumComplexes within Ostrich
ISEED integer Random seed used in optimization search. Enterany integer number. Default value (=1969) isassumed 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 theSCE-UA algorithm. Enter false or leave the fieldblank for default setting. Enter true for userspecified 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 shouldbe greater than or equal to 2. The defaultvalue is equal to (2 * number of optimizedparameters + 1).
---> this parameter is called m_PtsPerComplex within Ostrich
NPS integer Number of points in each sub-complex. NPSshould be greater than or equal to 2 and lessthan 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 complexbefore next shuffling. Default value is equalto NPG.
---> this parameter is called m_NumEvoSteps within Ostrich
MINGS integer Minimum number of complexes required foroptimization search, if the number of complexesis allowed to reduce as the optimization searchproceeds. 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 inthe 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 outthe best estimate of the global optimum at theend of each shuffling loop. Enter '1' to printout every point in the entire sample populationat the end of each shuffling loop. The defaultvalue is equal to 0. Enter 2 to ignore this variableand 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; //MAXNint m_Kstop; //KSTOPdouble m_Pcento; //PCENTOdouble m_Peps; //pepsdouble m_fSaved;int m_NumComplexes; //NGSint m_Seed; //ISEEDint m_UserConfig; //IDEFLTint m_PtsPerComplex; //NPGint m_PtsPerSubComplex; //NPSint m_NumEvoSteps; //NSPLint m_MinComplexes; //MINGSint m_UseInitPt; //INIFLGint m_OutputMode; //IPRINTint m_np; //number of parametersint m_CurIter;double * m_pParams; //PARAMSdouble * m_pLower; //LOWERdouble * m_pUpper; //UPPERbool 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 beforeoptimization is terminated. The purpose ofMAXN is to stop an optimization search beforetoo much computer time is expended. MAXNshould be set large enough so that optimizationis generally completed before MAXN trials areperformed. Recommended value is 10,000 (increase ordecrease as necessary).
---> this parameter is called m_Budget within Ostrich
KSTOP integer Number of shuffling loops in which the criterion must improve by the specifiedpercentage or else optimization will beterminated. Recommended value: 5.
---> this parameter is called m_Kstop within Ostrich
PCENTO double Percentage by which the criterion value mustchange 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 optimizationsearch. Minimum value is 1.Recommended value is between 2 and 20 dependingon the number of parameters to be optimized andon the degree of difficulty of the problem.
---> this parameter is called m_NumComplexes within Ostrich
ISEED integer Random seed used in optimization search. Enterany integer number. Default value (=1969) isassumed 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 theSCE-UA algorithm. Enter false or leave the fieldblank for default setting. Enter true for userspecified 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 shouldbe greater than or equal to 2. The defaultvalue is equal to (2 * number of optimizedparameters + 1).
---> this parameter is called m_PtsPerComplex within Ostrich
NPS integer Number of points in each sub-complex. NPSshould be greater than or equal to 2 and lessthan 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 complexbefore next shuffling. Default value is equalto NPG.
---> this parameter is called m_NumEvoSteps within Ostrich
MINGS integer Minimum number of complexes required foroptimization search, if the number of complexesis allowed to reduce as the optimization searchproceeds. 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 inthe 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 outthe best estimate of the global optimum at theend of each shuffling loop. Enter '1' to printout every point in the entire sample populationat the end of each shuffling loop. The defaultvalue is equal to 0. Enter 2 to ignore this variableand 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; //MAXNint m_Kstop; //KSTOPdouble m_Pcento; //PCENTOdouble m_Peps; //pepsdouble m_fSaved;int m_NumComplexes; //NGSint m_Seed; //ISEEDint m_UserConfig; //IDEFLTint m_PtsPerComplex; //NPGint m_PtsPerSubComplex; //NPSint m_NumEvoSteps; //NSPLint m_MinComplexes; //MINGSint m_UseInitPt; //INIFLGint m_OutputMode; //IPRINTint m_np; //number of parametersint m_CurIter;double * m_pParams; //PARAMSdouble * m_pLower; //LOWERdouble * m_pUpper; //UPPERbool 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 beforeoptimization is terminated. The purpose ofMAXN is to stop an optimization search beforetoo much computer time is expended. MAXNshould be set large enough so that optimizationis generally completed before MAXN trials areperformed. Recommended value is 10,000 (increase ordecrease as necessary).
---> this parameter is called m_Budget within Ostrich
KSTOP integer Number of shuffling loops in which the criterion must improve by the specifiedpercentage or else optimization will beterminated. Recommended value: 5.
---> this parameter is called m_Kstop within Ostrich
PCENTO double Percentage by which the criterion value mustchange 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 optimizationsearch. Minimum value is 1.Recommended value is between 2 and 20 dependingon the number of parameters to be optimized andon the degree of difficulty of the problem.
---> this parameter is called m_NumComplexes within Ostrich
ISEED integer Random seed used in optimization search. Enterany integer number. Default value (=1969) isassumed 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 theSCE-UA algorithm. Enter false or leave the fieldblank for default setting. Enter true for userspecified 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 shouldbe greater than or equal to 2. The defaultvalue is equal to (2 * number of optimizedparameters + 1).
---> this parameter is called m_PtsPerComplex within Ostrich
NPS integer Number of points in each sub-complex. NPSshould be greater than or equal to 2 and lessthan 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 complexbefore next shuffling. Default value is equalto NPG.
---> this parameter is called m_NumEvoSteps within Ostrich
MINGS integer Minimum number of complexes required foroptimization search, if the number of complexesis allowed to reduce as the optimization searchproceeds. 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 inthe 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 outthe best estimate of the global optimum at theend of each shuffling loop. Enter '1' to printout every point in the entire sample populationat the end of each shuffling loop. The defaultvalue is equal to 0. Enter 2 to ignore this variableand 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 filepFile = fopen(fileName, "a"); m_pStats->WriteStats(pFile);fclose(pFile);//write statistics of best parameter set to output filem_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 bannerWriteBanner(m_pModel, "gen best value ", "Pct. Complete");scemain(); //main SCE implemenation, converted from FORTRAN//place model at optimal prameter setpGroup->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 metricsWriteAlgMetrics(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 VARIABLESFILE * 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 PARAMETERSint 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 PARAMETERSif (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 BOUNDSint 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 VALUESif (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 VALIDif ((*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 POPULATIONint npt;npt = (*ngs) * (*npg); //npt = ngs * npgif (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 MESSAGESif (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);elsestrcpy(reduc, noflg);if (iniflg != 0)strcpy(initl, ysflg);elsestrcpy(initl, noflg);//PRINT SHUFFLED COMPLEX EVOLUTION OPTIMIZATION OPTIONSfprintf(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.1by QINGYUN DUAN (adpated to C++ by L. Shawn Matott)DEPARTMENT OF HYDROLOGY & WATER RESOURCESUNIVERSITY OF ARIZONA, TUCSON, AZ 85721(602) 621-9360, email: duan@hwr.arizona.eduWRITTEN IN OCTOBER 1990.REVISED IN AUGUST 1991REVISED IN APRIL 1992STATEMENT BY AUTHOR:--------------------This general purpose global optimization program is developed atthe Department of Hydrology & Water Resources of the Universityof Arizona. Further information regarding the SCE-UA method canbe obtained from Dr. Q. Duan, Dr. S. Sorooshian or Dr. V.K. Guptaat the address and phone number listed above. We request allusers of this program make proper reference to the paper entitled'Effective and Efficient Global Optimization for ConceptualRainfall-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 VARIABLESa(.) = initial parameter setbl(.) = lower bound on parametersbu(.) = upper bound on parametersnopt = number of parameters to be optimizedLIST OF SCE ALGORITHMIC CONTROL PARAMETERS:ngs = number of complexes in the initial populationnpg = number of points in each complexnpt = total number of points in initial population (npt=ngs*npg)nps = number of points in a sub-complexnspl = number of evolution steps allowed for each complex beforecomplex shufflingmings = minimum number of complexes required, if the number ofcomplexes is allowed to reduce as the optimization proceedsiseed = initial random seediniflg = flag on whether to include the initial point in population= 0, not included= 1, includediprint = 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 populationCONVERGENCE CHECK PARAMETERSmaxn = max no. of trials allowed before optimization is terminatedkstop = number of shuffling loops in which the criterion value mustchang by the given percentage before optimization is terminatedpcento = percentage by which the criterion value must change ingiven number of shuffling loopsipcnvg = 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 satisfiedLIST OF LOCAL VARIABLESx(.,.) = coordinates of points in the populationxf(.) = function values of x(.,.)xx(.) = coordinates of a single point in xcx(.,.) = coordinates of points in a complexcf(.) = function values of cx(.,.)s(.,.) = coordinates of points in the current simplexsf(.) = function values of s(.,.)bestx(.) = best point at current shuffling loopbestf = function value of bestx(.)worstx(.) = worst point at current shuffling loopworstf = function value of worstx(.)xnstd(.) = standard deviation of parameters in the populationgnrng = normalized geometric mean of parameter rangeslcs(.) = indices locating position of s(.,.) in x(.,.)bound(.) = bound on ith variable being optimizedngs1 = number of complexes in current populationngs2 = number of complexes in last populationiseed1 = current random seedcriter(.) = vector containing the best criterion values of the last10 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 ARRAYSdouble ** x, * xx, * bestx, * worstx, * xf;double ** s, * sf, ** cx, * cf;double * xnstd, * bound, criter[20], * unit;int * lcs; //allocate memorylcs = 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 VARIABLESint 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 INTEGERint iseed1;iseed1 = -abs(iseed);//COMPUTE THE TOTAL NUMBER OF POINTS IN INITIAL POPUALTIONint 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 OPTIMIZEDint 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 startif(m_pModel->CheckWarmStart() == true){WarmStart();m_pModel->GetParamGroupPtr()->ReadParams(a);}// handle parameter extractionif(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 budgetm_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 VALUEfprintf(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 SPACEif (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 oneif(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 VALUEsort(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 POPULATIONdouble gnrng;int ipcnvg;parstt(npt1,nopt,x,xnstd,bound,&gnrng,&ipcnvg);// PRINT THE RESULTS FOR THE INITIAL POPULATIONfprintf(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 COMPLEXESfor(igs = 1; igs <= ngs1; igs++){// ASSIGN POINTS INTO COMPLEXESint 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 DISTRIBUTIONif (nps == npg){int k;for(k = 0; k < nps; k++){lcs[k] = k;} // end dogoto label_85; } // end ifdouble 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 dolcs[k-1] = lpos; } // end do// ARRANGE THE SUB-COMPLEX IN ORDER OF INCEASING FUNCTION VALUEsort(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 dosf[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 COMPLEXfor(k = 1; k <= nps; k++) {for(j = 1; j <= nopt; j++) {cx[lcs[k-1]][j-1] = s[k-1][j-1]; } // end docf[lcs[k-1]] = sf[k-1]; } //end do// SORT THE POINTSsort(npg,nopt,cx,cf);//IF MAXIMUM NUMBER OF RUNS EXCEEDED, BREAK OUT OF THE LOOPif (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 doxf[k2-1] = cf[k1-1]; } // end doif (icall >= maxn) break; //END LOOP ON COMPLEXES} /* end for() */ // RE-SORT THE POINTSsort(npt1,nopt,x,xf);// RECORD THE BEST AND WORST POINTSfor(j = 0; j < nopt; j++) {m_pParams[j] = bestx[j] = x[0][j];worstx[j] = x[npt1-1][j]; } //end dom_Best = bestf = xf[0];worstf = xf[npt1-1];// TEST THE POPULATION FOR PARAMETER CONVERGENCEparstt(npt1,nopt,x,xnstd,bound,&gnrng,&ipcnvg);// PRINT THE RESULTS FOR CURRENT POPULATIONm_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 EXCEEDEDif (icall >= maxn) goto label_9000; // COMPUTE THE COUNT ON SUCCESSIVE LOOPS W/O FUNCTION IMPROVEMENTcriter[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 SPACEif (ipcnvg == 1) goto label_9200; //NONE OF THE STOPPING CRITERIA IS SATISFIED, CONTINUE SEARCH//CHECK FOR COMPLEX NUMBER REDUCTIONint 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 VALUEfprintf(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 VARIABLESconst double c1=0.8;const double c2=0.4;/* ----------------------------------------------------LIST OF LOCAL VARIABLESsb(.) = the best point of the simplexsw(.) = the worst point of the simplexw2(.) = the second worst point of the simplexfw = function value of the worst pointce(.) = the centroid of the simplex excluding wosnew(.) = new point generated from the simplexiviol = 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 CODEint n = nps;int m = nopt;double alpha = 1.0;double beta = 0.5;/* ---------------------------------------------------IDENTIFY THE WORST POINT wo OF THE SUB-COMPLEX sCOMPUTE THE CENTROID ce OF THE REMAINING POINTSCOMPUTE step, THE VECTOR BETWEEN wo AND ceIDENTIFY 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 doce[j] = ce[j]/(double)(n-1);} //end dodouble fw;fw = sf[n-1]; //COMPUTE THE NEW POINT snew//FIRST TRY A REFLECTION STEPfor(j = 0; j < m; j++) {snew[j] = ce[j] + alpha * (ce[j] - sw[j]);} //end do//CHECK IF snew SATISFIES ALL CONSTRAINTSint ibound;chkcst(nopt,snew,bl,bu,&ibound); /* ------------------------------------------------------------------snew IS OUTSIDE THE BOUND,CHOOSE A POINT AT RANDOM WITHIN FEASIBLE REGION ACCORDING TOA NORMAL DISTRIBUTION WITH BEST POINT OF THE SUB-COMPLEXAS 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 snewdouble 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 RETURNif (fnew <= fw) goto label_2000;if (*icall >= maxn) goto label_3000;//fnew IS GREATER THAN fw, SO TRY A CONTRACTION STEPfor(j = 0; j < m; j++){snew[j] = ce[j] - beta * (ce[j] - sw[j]);} //end do//COMPUTE THE FUNCTION VALUE OF THE CONTRACTED POINTeb = (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 RETURNif (fnew <= fw) goto label_2000; if (*icall >= maxn) goto label_3000; /* ---------------------------------------------------------------------IF BOTH REFLECTION AND CONTRACTION FAIL, CHOOSE ANOTHER POINTACCORDING TO A NORMAL DISTRIBUTION WITH BEST POINT OF THE SUB-COMPLEXAS 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 POINTeb = (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 dosf[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 constraintschkcst(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 VALUESdouble 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 doxmean[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 ifdelete [] 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 dobf[k2-1] = af[k1-1];} //end do} //end doint j, jj;for(j = 0; j < npt; j++) {jj = j % npg;for(i = 0; i < n; i++){a[j][i] = b[jj][i]; } //end doaf[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 costsfor(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 pthif(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 pthif(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 allconstraints.ibound - violation indicator= -1 initial value= 0 no violation= 1 violationnopt = number of optimizing variablesii = 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 violatedfor(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 defaultsm_np = m_pModel->GetParamGroupPtr()->GetNumParams();m_Budget = 10000; //MAXNm_Kstop = 5; //KSTOPm_Pcento = 0.01; //PCENTOm_Peps = 1.0E-3; //pepsm_NumComplexes = (int)(sqrt((double)m_np)); //NGSm_Seed = 1969; //ISEEDm_UserConfig = 1; //IDEFLTm_PtsPerComplex = 2*m_np + 1; //NPGm_PtsPerSubComplex = m_np+1; //NPSm_NumEvoSteps = m_PtsPerComplex; //NSPLm_MinComplexes = m_NumComplexes; //MINGSm_UseInitPt = 1; //INIFLGm_OutputMode = 2; //IPRINTm_bUseInitPt = false; //INIFLG//allocate initial parameter configurationParameterGroup * 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 configurationpFile = 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 presentif(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 fileint 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() */
【算法】02 SCE-UA简介及源代码相关推荐
- C#，图像二值化（20）——全局阈值的耶恩算法（Yen Thresholding）及源代码
1 全局阈值的耶恩算法(Yen Throsholding) 常见阈值算法 1.1黄算法 HuangThresholdImageFilter使用Shannon的熵函数实现Huang的模糊阈值[1].模糊 ...
- 零起点学算法02——输出简单的句子
零起点学算法02--输出简单的句子 Description 会输出Hello World!了,那换个句子也会吧? Input 没有输入 Output 现在要求你输出下面红色的字 Nice to me ...
- 弗洛伊德算法（Floyd）简介
弗洛伊德算法(Floyd)简介 Floyd算法适用于解决求最短路径问题,相比于Digkstra算法思路更加简单,更容易理解,但是效率会明显低很多,可以作为初步学习的一种方法. 目录 弗洛伊德算法(Fl ...
- Interview之NLP：人工智能领域求职岗位—自然语言处理NLP算法工程师职位的简介、薪资介绍、知识结构之详细攻略
Interview之NLP:人工智能领域求职岗位-自然语言处理NLP算法工程师职位的简介.薪资介绍.知识结构之详细攻略 目录 自然语言处理NLP算法工程师的职位简介 1.资讯指数 2.各大公司的具体职 ...
- Interview之CV：人工智能领域求职岗位—计算机视觉算法工程师的职位简介、薪资介绍、知识结构之详细攻略
Interview之CV:人工智能领域求职岗位-计算机视觉算法工程师的职位简介.薪资介绍.知识结构之详细攻略 目录 计算机视觉算法工程师的职位简介 资讯指数 1.各大互联网巨头的薪资介绍 2.各大公司 ...
- 手把手的K-means聚类算法教程（含简介及教育数据应用实例 Python实现）
手把手的K-means聚类算法教程(含简介及教育数据应用实例 Python实现) 1. K-MEANS的基本原理 2. 数据预处理 2.1 数据读取:from Excel 2.2 数据预处理:标准化Z ...
- 结构与算法(02)：队列和栈结构
本文源码:GitHub·点这里 || GitEE·点这里 一.队列结构 1.基础概念 队列是一种特殊的线性表,特殊之处在于它只允许在表的前端(front)进行删除操作,而在表的后端(rear)进行插入 ...
- C#，图像二值化（16）——全局阈值的力矩保持算法（Moment-proserving Thresholding）及其源代码
1.力矩保持法 提出了一种基于矩保持原理的自动阈值选择方法.以这样的方式确定地计算阈值,即在输出画面中保留输入画面的时刻.实验结果表明,该方法可以将给定的图像阈值化为有意义的灰度级.该方法描述了全局阈 ...
- C#，图像二值化（12）——基于谷底最小值的全局阈值算法（Valley-Minium Thresholding）与源代码
1.基于谷底最小值的阈值 这个方法适用于通过有限的迭代次数,平滑后能得到双峰的图像,让双峰的谷底成为阈值.当执行完基于谷底最小值的阈值操作的时候,改变了直方图信息,使之成为处理过后的直方图信息,这时候 ...
- C#，图像二值化（06）——全局阈值的大津算法（OTSU Thresholding）及其源代码
1.大津算法OTSU ALGORITHM OTSU算法效果很一般. 最大类间方差法是1979年由日本学者大津(Nobuyuki Otsu)提出的,是一种自适应阈值确定的方法,又叫大津法,简称OTSU, ...
最新文章
- C++中Reference与指针（Pointer）的使用对比
- svn冲突的解决办法
- Erlang error？
- scikit-learn kmeans++
- Entity Framework 学习结束语
- 实战SSM_O2O商铺_06logback配置与使用
- 使用ilmerge实现.net程序静态链接
- java判断表是否存在_java怎么判断表是否存在？
- clickhouse原理解析与开发实战 pdf_重识SSM，“超高频面试点+源码解析+实战PDF”，一次性干掉全拿走...
- python性能解决_Python性能优化的20条建议
- King of the Ether
- 一种抑制undershoot/overshoot锐化算法介绍
- js 20160810
- 超简单的Tomcat安装过程
- 本周新出开源计算机视觉代码汇总（含图像超分辨、视频目标分割、行人重识别、点云识别等）...
- readonly和disabled的区别
- Linux下运行java DES AES加解密
- 想学PHP来兄弟连是正确的选择 初识兄弟连三周
- 网易云爬取歌词进行歌词词云可视化
- php字游戏源码,php文字游戏寻仙纪.zip