C               OPTIMISATION ROUTINE UTILISING

C                    A GENETIC 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 This 'main' program initializes the GA program

C

      IMPLICIT REAL (A-H,O-Z)

C

      INTEGER J,K,N,SL,NC,IT,AL,RST

      INTEGER SOL(50,30,2)

      REAL PC,PM

C

C Read in data

C

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

C

      READ(5,*) NSEED

C

C Initialize random number generator

C

      CALL SEEDRN(NSEED)

C

C Population size:

C

      READ(5,*) N

C

C String length:

C

      READ(5,*) SL

C

C Number of chromosomes:

C

      READ(5,*) NC

C

C Number of generations

C

      READ(5,*) IT

C

C Alphabet size

C

      READ(5,*) AL

C

C The crossover probability

C

      READ(5,*) PC

C

C The mutation probability

C

      READ(5,*) PM

C

C Is this a restarted search?

C RST = 1 if so, RST = 0 if not

C

      READ(5,*) RST

C

C If not...

C Generate initial population randomly

C

      IF(RST .EQ. 0) THEN

         CALL INIPOP(N,SL,NC,AL,SOL)

      ELSE

C

C Otherwise read in population

C SOL(i,j) holds the allele value

C

         DO 1 K = 1,N

            DO 1 KK = 1,NC

               READ(5,*) (SOL(K,J,KK),J=1,SL)

 1       CONTINUE

      ENDIF

      CLOSE(UNIT=5)

C

C Now call GA program

C

      CALL GA(N,SL,NC,IT,AL,PC,PM,SOL)

C

      END

C

C

      SUBROUTINE GA(N,SL,NC,IT,AL,PC,PM,SOL)

      INTEGER N,SL,NC,IT,AL,I,J,K

      INTEGER SOL(50,30,2),NEWS(50,30,2)

      INTEGER PRS(25,2),RK(50),BSOL

      REAL PC,PM,FMAX,FAVE

      REAL FI(50),SP(50)

      COMMON / BSOLN / BSOL(30,2),FBEST]

C

C Set the initial best objective function to a suitably large value

C 

      FBEST = -1.E32

C  

C The heading for generation by generation monitoring

C

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

      WRITE(6,992) ' GEN     FMAX      FAVE '

C

C The main loop of the GA itself

C

      DO 1 I = 1,IT

C

C Assess the current population - find the fitness of each solution

C

         CALL MODEL(I,N,SL,NC,SOL,FI)

C

C Rank the solutions according to fitness - in order to use Baker's

C survival probability criterion - and find the maximum and average

C values

C

         CALL ORDER(N,FI,RK,FMAX,FAVE,IMAX)

C

C Generation by generation monitoring

C

         WRITE(6,990) I,FMAX,FAVE

C

C Find the survival probability for each solution using Baker's

C ranking rule and put it in SP

C

         CALL SURVIV(N,RK,SP)

C

C If on the last iteration, do not perturb the assessed population

C

         IF(I .EQ. IT) GOTO 1

C

C Select the new population using 'stochastic remainder selection 

C without replacement'

C

         CALL SELECT(N,SL,NC,SOL,SP,NEWS,FMAX,IMAX)

C

C Randomly pair off the new population prior to crossover

C

         CALL PAIRS(N,PRS)

C

C Perform crossover between these pairs

C

         CALL CROSS(N,SL,NC,PC,NEWS,PRS)

C

C Perform mutation on new solutions

C

         CALL MUTATE(N,SL,NC,AL,PM,NEWS,SOL)

C

 1    CONTINUE

C

C End of run output

C

      WRITE(6,*)

      WRITE(6,991) 'LAST GENERATION:'

      WRITE(6,*)

      WRITE(6,992) ' SOL     FI     RK    SP'

      DO 2 K = 1,N

         WRITE(6,993) K,FI(K),RK(K),SP(K)

 2    CONTINUE

      CLOSE(UNIT=6)

C

C Output the data necessary for a restart

C

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

      NSEED = INT(1000*RANGEN())

      WRITE(7,996) NSEED

      WRITE(7,996) N

      WRITE(7,996) SL

      WRITE(7,996) NC

      WRITE(7,996) IT

      WRITE(7,996) AL

      WRITE(7,997) PC

      WRITE(7,997) PM

      WRITE(7,998) '   1'

      DO 3 K = 1,N

         DO 3 KK = 1,NC

            WRITE(6,994) (SOL(K,J,KK),J=1,SL)

 3    CONTINUE

      CLOSE(UNIT=7)

 990  FORMAT(I4,2F10.4)

 991  FORMAT(A17)

 992  FORMAT(A24)

 993  FORMAT(I4,F10.4,I4,F8.4)

 994  FORMAT(40I1)

 995  FORMAT(A18)

 996  FORMAT(I4)

 997  FORMAT(F6.3)

 998  FORMAT(A4)

C

      RETURN

      END

C

C This subroutine generates a random initial population

C

      SUBROUTINE INIPOP(N,SL,NC,AL,SOL)

      INTEGER N,SL,NC,AL,I,J,K

      INTEGER SOL(50,30,2)

      DO 1 I = 1,N

         DO 1 J = 1,SL

            DO 1 K = 1,NC

               SOL(I,J,K) = AL*RANGEN()

 1    CONTINUE

      RETURN

      END

C

C This subroutine ranks the population according to fitness

C 1st = best, nth = worst

C To avoid ties between solutions with the same fitness that with

C the lower index is deemed to be better

C It also finds the maximum and average fitness

C

      SUBROUTINE ORDER(N,FI,RK,FMAX,FAVE,IMAX)

      INTEGER N,I,J,R

      INTEGER RK(50)

      REAL FI(50),FMAX,FAVE

      FMAX = -1.E32

      FAVE = 0.0

      DO 1 I = 1,N

         IF(FI(I).GT.FMAX) THEN

            FMAX = FI(I)

            IMAX = I

         ENDIF 

         FAVE = FAVE + FI(I)

         R = 1

         DO 2 J = 1,N

            IF(I .LT. J) THEN

               IF(FI(I) .LT. FI(J)) R = R + 1

            ENDIF

            IF(I .GT. J) THEN

               IF(FI(I) .LE. FI(J)) R = R + 1

            ENDIF

 2       CONTINUE

         RK(I) = R

 1    CONTINUE

      FAVE = FAVE/N

      RETURN

      END

C

C This subroutine calculates the survival probabilities according

C to Baker's ranking criterion

C

      SUBROUTINE SURVIV(N,RK,SP)

      INTEGER N,I

      INTEGER RK(50)

      REAL SP(50)

      DO 1 I = 1,N

         SP(I) = (2.000001*(N-RK(I)))/(N*(N-1.0))

 1    CONTINUE

      RETURN

      END

C

C This subroutine selects the new population according by

C `stochastic remainder selection without replacement'

C

      SUBROUTINE SELECT(N,SL,NC,SOL,SP,NEWS,FMAX,IMAX)

      INTEGER N,SL,NC,I,J,K,CP,NA,BSOL

      INTEGER SOL(50,30,2),NEWS(50,30,2)

      REAL SP(50),E(50),R(50)

      REAL SEL

      COMMON / BSOLN / BSOL(30,2),FBEST

C

C If the best solution in the current generation is better than 

C the recorded best, update the record

C

      IF(FMAX.GE.FBEST) THEN

         FBEST = FMAX

         DO 22 J = 1,SL

            DO 22 K = 1,NC

               BSOL(J,K) = SOL(IMAX,J,K)

 22      CONTINUE

      ENDIF

C

C Insert a copy of the best solution in the new population

C

      CP=1

      DO 23 J=1,SL

         DO 23 K = 1,NC

            NEWS(CP,J,K)=BSOL(J,K)

 23   CONTINUE

C

C Calculate the expected number of copies of each solution

C and make those that are assured

C

      DO 1 I = 1,N

         E(I) = N*SP(I)

         NA = INT(E(I))

C

C Make one less copy of the best solution

C

         IF(I.EQ.IMAX) NA = NA-1

         R(I) = E(I) - NA

         DO 2 J = 1,NA

            CP = CP + 1

            DO 3 K = 1,SL

               DO 3 KK = 1,NC

                  NEWS(CP,K,KK) = SOL(I,K,KK)

 3          CONTINUE

 2       CONTINUE

 1    CONTINUE

C

C Now make the rest by probabilistic selection

C

 4    IF(CP .EQ. N) GOTO 10

      DO 5 I = 1,N

         SEL = RANGEN()

         IF(SEL .LT. E(I)) THEN

            CP = CP + 1

            DO 6 K = 1,SL

               DO 6 KK = 1,NC

               NEWS(CP,K,KK) = SOL(I,K,KK)

 6          CONTINUE

         ENDIF

         IF(CP .EQ. N) GOTO 10

 5    CONTINUE

      IF(CP .LT. N) GOTO 4

 10   RETURN

      END

C

C This subroutine randomly pairs off the population

C

      SUBROUTINE PAIRS(N,PRS)

      INTEGER I,J,K,N,X,Y

      INTEGER PRS(25,2),TKN(50)

C

C The array TKN holds the solutions not yet taken

C

      DO 1 I = 1,N

         TKN(I) = I

 1    CONTINUE

      DO 2 I = 1,N

         J = 1+ (N+1-I)*RANGEN()

         X = 1 + (I-1)/2

         Y = 2 - MOD(I,2)

         PRS(X,Y) = TKN(J)

         DO 3 K = J,N-1

            TKN(K) = TKN(K+1)

 3       CONTINUE

 2    CONTINUE

      RETURN

      END

C

C This subroutine performs two point crossover between each

C pair of solutions with probability PC

C

      SUBROUTINE CROSS(N,SL,NC,PC,NEWS,PRS)

      INTEGER I,J,F,M,N,SL,NC,X,Y

      INTEGER NEWS(50,30,2),PRS(25,2),BUF

      REAL PC,SEL

      DO 1 I = 1,N/2

         SEL = RANGEN()

C

C If crossover is approved, randomly select two coordinates X & Y

C

         IF(SEL .LT. PC) THEN

C

C For each chromosome...

C 

           DO 11 KK = 1,NC

              F = PRS(I,1)

              M = PRS(I,2)

              X = 1+ SL*RANGEN()

              Y = 1+ SL*RANGEN()

              IF(Y .GE. X) THEN

                 DO 2 J = X,Y

                    BUF = NEWS(F,J,KK)

                    NEWS(F,J,KK) = NEWS(M,J,KK)

                    NEWS(M,J,KK) = BUF

 2             CONTINUE

            ENDIF

            IF(Y .LT. X) THEN

               DO 4 J = X,SL

                  BUF = NEWS(F,J,KK)

                  NEWS(F,J,KK) = NEWS(M,J,KK)

                  NEWS(M,J,KK) = BUF

 4             CONTINUE

               DO 6 J = 1,Y

                  BUF = NEWS(F,J,KK)

                  NEWS(F,J,KK) = NEWS(M,J,KK)

                  NEWS(M,J,KK) = BUF

 6             CONTINUE

            ENDIF

 11        CONTINUE

         ENDIF

 1    CONTINUE

      RETURN

      END

C

C This subroutine performs mutation on solutions with probability

C PM and returns solutions in array SOL ready for evaluation

C

      SUBROUTINE MUTATE(N,SL,NC,AL,PM,NEWS,SOL)

      INTEGER N,SL,NC,AL,I,J,X

      INTEGER NEWS(50,30,2),SOL(50,30,2)

      REAL PM,SEL

      DO 1 I = 1,N

         DO 1 J = 1,SL

            DO 1 KK = 1,NC 

               SEL = RANGEN()

               IF(SEL .LT. PM) THEN

 3                X = AL*RANGEN()

                  IF(X .EQ. NEWS(I,J,KK)) GOTO 3

                  NEWS(I,J,KK) = X

               ENDIF

            SOL(I,J,KK) = NEWS(I,J,KK)

 1    CONTINUE

      RETURN

      END

C

      SUBROUTINE MODEL(IGEN,N,SL,NC,SOL,FI)

C

C Rosenbrock's function for 2 variables

C

      INTEGER N,SL,NC

      INTEGER SOL(50,30,2)

      REAL FI(50),X(50,2),RX(2)

      WRITE(66,*) IGEN

C

C Convert solutions into real numbers between 0 & 1

C

      DO 1 I = 1,N

         DO 1 KK = 1,NC

           X(I,KK) = 0.0

           DO 2 J = 1,SL

              IF(SOL(I,J,KK).EQ.1)

     &           X(I,KK) = X(I,KK) + (2.0**(J-1))/(2.0**SL)

 2        CONTINUE

 1    CONTINUE

C

C Convert these numbers into rescaled reals and evaluate the

C objective function

C These scalings are 'hard-wired': RX1 varies from -1.5 to 1.5

C                                  RX2 varies from -1.0 to 2.0 

C

      DO 3 I = 1,N

         RX(1) = 3*X(I,1) - 1.5

         RX(2) = 3*X(I,2) - 1.0

         FI(I) = (1.0-RX(1))**2 + 100.0*(RX(2)-RX(1)*RX(1))**2

C

C As the GA is a maximization routine negate the objective value

C

         FI(I) = -1*FI(I)

 3    CONTINUE

      RETURN

      END
