C            OPTIMISATION ROUTINE UTILISING THE

C        METROPOLIS (SIMULATED ANNEALING) ALGORITHM

C                A STOCHASTIC SEARCH PROGRAM

C        UNCONSTRAINED CONTINUOUS VARIABLES VERSION

C

C The user must supply two machine-dependent subroutines/functions:

C  SEEDRN(nseed) a subroutine to seed the random number generator

C  RANGEN() a function returning random numbers between 0 and 1

C

C Declare all variables used

C

      IMPLICIT REAL (A-H,O-Z)

C

      INTEGER NV,I,J,INC,REP,RED,ARD,NT,INX

      INTEGER INCS,INXS,REDS,ARDS,ITS

      INTEGER LTRAN,MACCS

C

      REAL OF,OFO,DOF,OFMIN,SDOF,DFAVE

      REAL STAVE,T,T1,PLIM,P,RAT,ALFA

C

      REAL RX(10),RXL(10),RXU(10),X(10)

      REAL XO(10),DX(10),XMIN(10),STMX(10)

      REAL STMXL(10),STMXU(10)

C

C Read in data

C

      OPEN(UNIT=5,FILE='SA.DAT')

C

C Random number generator seed

C

      READ(5,*) NSEED

C

C Initialize random number generator

C

      CALL SEEDRN(NSEED)

C

C Number of control variables

C

      READ(5,*) NV

C

C Initial values of control variables - real scaling

C Lower bounds on control variables - real scaling

C Upper bounds on control variables - real scaling

C

      DO 1 I = 1,NV

         READ(5,*) RX(I),RXL(I),RXU(I)

 1    CONTINUE

C

C Scale control variables so that they vary between -1 and +1

C

      DO 2 I = 1,NV

         X(I) = (2*RX(I)-RXU(I)-RXL(I))/(RXU(I)-RXL(I))

 2    CONTINUE

C

C Initial estimates/values of appropriate step sizes

C Lower bounds on step sizes

C Upper bounds on step sizes

C

      STAVE = 0.0

      DO 3 I = 1,NV

         READ(5,*) STMX(I),STMXL(I),STMXU(I)

         STAVE = STAVE + STMX(I)/NV

 3    CONTINUE

C

C Control parameters for algorithm

C

C Maximum number of trials to be performed

C

      READ(5,*) NT

C

C Maximum number of trials at a given temperature

C

      READ(5,*) LTRAN

C

C Maximum number of acceptances at a given temperature

C

      READ(5,*) MACCS

C

C Initial system temperature -

C if set to zero T1 iteration is required

C

      READ(5,*) T1

C

C Cooling rate control parameter

C

      READ(5,*) ALFA

      CLOSE(UNIT=5)



C

C Initialize performance monitoring counters

C

C Total number of objective function increases generated

C

      INX = 0

      INXS = 0

C

C Number of objective function increases accepted

C

      INC = 0

      INCS = 0

C

C Number of objective function reductions generated

C

      RED = 0

      REDS = 0

C

C Number of absolute objective function reductions found

C

      ARD = 0

      ARDS = 0

C

C Initialize performance controlling counters

C

C Number of trials since successful one -

C increase or decrease of objective function

C Number of consecutive 'complete failures'

C

      REP = 0

C

C Number of trials at current temperature

C

      ITS = 0

C

C Print heading for output

C

      OPEN(UNIT=6,FILE='MON.DAT')

      WRITE(6,*) '  N       OFO     OFMIN     INC   RED   ARD  ',

     &           ' STAVE      T         RAT'

C

C Call subroutine MODEL to find initial value of objective function

C

      CALL MODEL(NV,X,RXL,RXU,OF)

C

C Store starting point as best solution so far

C

      DO 4 J = 1,NV

         XO(J) = X(J)

         XMIN(J) = X(J)

 4    CONTINUE

       OFO = OF

       OFMIN = OF

C

C If initial temperature estimation is required, conduct a 100 trial

C survey

C

      IF(T1 .EQ. 0.0) THEN

         SDOF = 0.0

         DO 44 I = 1,100

C

C Generate random new position ...

C

            DO 45 J = 1,NV

               XXX = 2*RANGEN()-1.0

               X(J) = XO(J) + STMX(J)*XXX

C

C Restrict X(J) to lie between -1 and +1

C

               IF(X(J) .GT. 1.0) X(J) = 1.0

               IF(X(J) .LT. -1.0) X(J) = -1.0

 45         CONTINUE

C

C Evaluate objective function at new point

C

            CALL MODEL(NV,X,RXL,RXU,OF)

C

C If the objective function is the smallest so far

C

            IF(OF .LT. OFMIN) THEN

C

C Then update stores

C

               OFMIN = OF

               DO 46 J = 1,NV

                  XMIN(J) = X(J)

 46            CONTINUE

                ARD = ARD + 1

            ENDIF

C

C If the objective is increased, update statistics

C

            IF(OF .GT. OFO) THEN

               INX = INX + 1

               DOF = OF - OFO

               SDOF = SDOF + DOF

            ENDIF

C

C In this initial survey, accept all changes

C

            DO 47 J = 1,NV

               XO(J) = X(J)

 47         CONTINUE

             OFO = OF

 44      CONTINUE

C

C At end of survey calculate appropriate initial temperature

C (To give initial acceptance ratio of 0.8)

C (ln(1.25) = 0.22314)

C

         DFAVE = SDOF/INX

         T1 = 0.22314*DFAVE

         INX = 0

      ENDIF

C

C Start of algorithm itself

C

C Set initial temperature

C

      T = T1

C

C Randomly search through variable space seeking improvement

C

      DO 5 I = 1,NT

         ITS = ITS + 1

C

C Generate random new position ...

C

C Use self-tuning step size for each variable

C

         DO 8 J = 1,NV

            XXX = 2*RANGEN()-1.0

            X(J) = XO(J) + STMX(J)*XXX

C

C Restrict X(J) to lie between -1 and +1

C

C If X(J) is at either limit then DX(J) will probably be greater

C than the actual change effected

C

            IF(X(J) .GT. 1.0) X(J) = 1.0

            IF(X(J) .LT. -1.0) X(J) = -1.0

            DX(J) = X(J) - XO(J)

 8       CONTINUE

C

C Evaluate objective function at new point

C

         CALL MODEL(NV,X,RXL,RXU,OF)

C

C If the objective function is the smallest so far

C

         IF(OF .LT. OFMIN) THEN

C

C Then update stores

C

            OFMIN = OF

            DO 9 J = 1,NV

               XMIN(J) = X(J)

 9          CONTINUE

             ARD = ARD + 1

         ENDIF

C

C If objective function is reduced from current value,

C accept new position

C

         IF(OF .LT. OFO) THEN

             OFO = OF

            DO 10 J = 1,NV

               XO(J) = X(J)

 10         CONTINUE

             RED = RED + 1

             REP = 0

         ELSE

C

C If OF is increased, accept new position with probability

C dependent on size of increase

C

C Change in objective function

C

            DOF = OF - OFO

            INX = INX + 1

C

C Apply the Metropolis test

C

            PLIM = EXP(-DOF/(STAVE*T))

            P = RANGEN()

            IF(P .LT. PLIM) THEN

                OFO = OF

               DO 11 J = 1,NV

                  XO(J) = X(J)

 11            CONTINUE

                INC = INC + 1

                REP = 0

            ELSE

               REP = REP + 1

            ENDIF

         ENDIF

C

C If trial has been a success update individual STMX settings

C There are upper and lower limits on the step sizes

C

         IF(REP .EQ. 0) THEN

             STAVE = 0.0

            DO 12 J = 1,NV

               STMX(J) = 0.9*STMX(J) + 0.21*ABS(DX(J))

               IF(STMX(J) .GT. STMXU(J)) STMX(J) = STMXU(J)

               IF(STMX(J) .LT. STMXL(J)) STMX(J) = STMXL(J)

               STAVE = STAVE + STMX(J)/NV

 12         CONTINUE

         ENDIF

C

C Decrease temperature after LTRAN trials or after MACCS acceptances

C

         IF((MOD(ITS,LTRAN) .EQ. 0) .OR. 

     &      ((INC+RED) .GT. MACCS)) THEN

            RAT = 1.0*INC/INX

            INCS = INCS + INC

            REDS = REDS + RED

            ITS = 0

            INC = 0

            INX = 0

            RED = 0

C

C Check for convergence - no absolute reduction at this temperature

C and acceptance ratio less than 0.1

C

            IF((ARD .EQ. ARDS) .AND. (RAT .LT. 0.10)) GOTO 665

            ARDS = ARD

C

C Output monitoring information

C

            WRITE(6,1112) I,OFO,OFMIN,INCS,REDS,ARD,STAVE,T,RAT

 1112       FORMAT(I6,2F10.4,3I6,F9.4,F11.4,F9.4)

C

C Temperature is actually changed here

C

            T = ALFA*T

         ENDIF

C

C End of one iteration

C

 5    CONTINUE

C

C Final output information

C

 665  WRITE(6,1112) I,OFO,OFMIN,INCS,REDS,ARD,STAVE,T,RAT

C

C Unscale control variables at best solution

C

      DO 13 J = 1,NV

         RX(J) = (XMIN(J)*(RXU(J)-RXL(J))+RXU(J)+RXL(J))/2

 13   CONTINUE

C

C Output best solution found

C

      WRITE(6,*)

      WRITE(6,*) 'BEST SOLUTION FOUND:'

      WRITE(6,1114) 'Objective Function:',OFMIN

      DO 14 J = 1,NV

         WRITE(6,1115) 'Variable ',J,':',RX(J)

 14   CONTINUE

      CLOSE(UNIT=6)

 1114 FORMAT(A20,F10.4)

 1115 FORMAT(A10,I2,A1,F9.4)

C

C Output the data necessary for a restart

C

      OPEN(UNIT=7,FILE='RST.DAT')

      NSEED = INT(1000*RANGEN())

      WRITE(7,1118) NSEED

      WRITE(7,1118) NV

      DO 16 J = 1,NV

         WRITE(7,1119) RX(J),RXL(J),RXU(J)

 16   CONTINUE

      DO 17 J = 1,NV

         WRITE(7,1119) STMX(J),STMXL(J),STMXU(J)

 17   CONTINUE

      WRITE(7,1118) NT

      WRITE(7,1118) LTRAN

      WRITE(7,1118) MACCS

      WRITE(7,1117) T1

      WRITE(7,1117) ALFA

      CLOSE(UNIT=7)

 1117 FORMAT(D13.5)

 1118 FORMAT(I6)

 1119 FORMAT(3D13.5)

C

 666  STOP

      END

C

      SUBROUTINE MODEL(NV,X,RXL,RXU,OF)

C

C Rosenbrock's function for NV variables

C

      INTEGER I,NV

      REAL X(10),RX(10),RXL(10),RXU(10),OF

C

C Scale the control variables from the range (-1,1) to (RXL,RXU)

C

      DO 1 I = 1,NV

        RX(I) = (X(I)*(RXU(I)-RXL(I))+RXU(I)+RXL(I))/2

 1    CONTINUE

C

C Calculate Rosenbrock's function

C

      OF = 0.0

      DO 2 I = 1,NV,2

         OF = OF + (1.0-RX(I))**2 

     &        + 100.0*(RX(I+1)-RX(I)*RX(I))**2

 2    CONTINUE

      RETURN

      END