!******************************************************************************** !* Genetic Algorithm(GA) program to minimize a function of two variables. !* !* written by Dr. Katate Masatsuka (info[at]cfdbooks.com), !* !* the author of useful CFD books, "I do like CFD" (http://www.cfdbooks.com). !* !* This is Version 1 (2011). !* !* This program was originally written in 1999 when I was interested in !* aerodynamic optimization problems. !* !* This file may be updated in future. !* !* Katate Masatsuka, January 2011. http://www.cfdbooks.com !******************************************************************************** !******************************************************************************** !* Genetic Algorithm(GA) for a function Minimization. !* (with "Npar" Continuous Parameters) !* !* This program finds the global minimizer of fnc(x1,...,xn). n=Npar. !* !* 1. Produce initial population randomly. ("Nipop" Chromosomes) !* 2. Select good chromosomes for GA. ("Npop " Chromosomes) !* 3.>>>>> Genetic Algorithm Loop >>>>> (Maximum # of loop = Imax) !* (1). Evaluate the cost for "Npop" Chromosomes. !* (2). Sort them in the decreasing order. !* (The upper "Nkeep" chromosomes are therefore considered as excellent) !* (3). Check for convergence. If converged, EXIT the loop. !* (4). Mating among "Nkeep" chromosomes to produce "Npop-Nkeep" offsprings. !* ("Npop-Nkeep" chromosomes will be replaced by the new offsprings.) !* (5). Mutation !* <<<<< Return to the top of Genetic Algorithm Loop <<<<< !* 4. End of the program. !* !* NB: This program uses very basic algorithms. It is intended only for !* an educational use. !* !* Katate Masatsuka, January 2011. http://www.cfdbooks.com !******************************************************************************** program genetic_algorithm implicit none !Numeric parameters integer , parameter :: sp = kind(1.0) integer , parameter :: dp = selected_real_kind(2*precision(1.0_sp)) real(dp), parameter :: zero=0.0_dp, one=1.0_dp, ten = 10.0_dp !Local parameters integer, parameter :: npar=2 !Parameters integer, parameter :: npop=20 !Tolal population integer, parameter :: nkeep=10 !Population to keep: the rest (npop-nkeep) ! is reserved for offsprings. !Local variables real(dp) :: chrom(npop,npar+1) !chrom(i,1:3)=[parameter01,parameter02,cost]_i real(dp) :: pmx(2,npar) !Min/max of the parameters real(dp) :: minc, meanc, stdc, tol_std !Min, mean, STD, tolerance integer :: index(npop) !Ordered index integer :: i, k, nparent1, nparent2 !***************************************************************************** !* Tolerance on standard deviation. tol_std = 1.0e-13_dp !***************************************************************************** !* Parameter ranges pmx(1,1) = zero pmx(2,1) = ten pmx(1,2) = zero pmx(2,2) = ten !***************************************************************************** !* Produce initial population: randomly produce "nipop" chromosomes call initial(chrom,npop,npar,pmx) do i = 1, npop index(i) = i end do call heapsorti(chrom,index,npop,-1,npar) !***************************************************************************** !* Genetic algorithm starts here. !***************************************************************************** ga_loop : do i = 1, 50000 ! Max number of revolution is 50000 (should be enough). !------------------------------------------------------------------ ! Evaluate the cost and sort the generation i for i>1. if (i/=1) then do k=1,npop chrom(k,npar+1) = cost_func(chrom,k,npar) end do call heapsorti(chrom,index,npop,-1,npar) endif !------------------------------------------------------------------ ! Evaluate the cost statistics over "nkeep" chromosomes: minc, meanc, stndc meanc = sum(chrom(index(1:nkeep),npar+1))/dble(nkeep) stdc = stndd(chrom,index,npar,nkeep,meanc) minc = minval(chrom(:,npar+1)) write(*,'(a14,i4,3(2x,a4,es12.2))') " Generation i=",i, & "min=",minc,"max=",meanc,"std=",stdc !----------------------------------------------------------------- ! Check for convergence (check over "nkeep" chromosomes only) if (stdc < tol_std) exit ga_loop !------------------------------------------------------------------ ! Pair chromosomes (from "npop" chrom) and produce "npop-nkeep" offsprings. do k = npop-nkeep+1, npop, 2 nparent1 = tournament(chrom,index,nkeep,npar) nparent2 = tournament(chrom,index,nkeep,npar) call crossover(chrom,index,k,nparent1,nparent2,npar) end do !------------------------------------------------------------------ ! Mutations call mutation(chrom,index,npop,nkeep,npar,pmx) end do ga_loop write(*,*) write(*,'(a37,es10.2)') ">>> GA converged with STD tolerance = ",tol_std write(*,'(a17,i5)') " n_population = ", npop write(*,'(a17,i5)') " n_keep = ", nkeep write(*,'(a17,i5)') " n_offspring = ", npop-nkeep !***************************************************************************** !* genetic algorithm ends here. !***************************************************************************** !***************************************************************************** !* display the final generation !***************************************************************************** write(*,*) write(*,*) "******************** Final Generation ***********************" do k=1,nkeep write(*,'(i6,3es16.3)') k, chrom(index(k),1:npar+1) end do write(*,*) "*************************************************************" !***************************************************************************** contains !***************************************************************************** !* Cost function to be minimized: A function of (x(1),...,x(npar)) !* !* Below, a simple function is set (the minimum is zero at x(1)=x(2)=0). !* Modify this function to try other functions. !**************************************************************************** function cost_func(chrom,i,npar) result(cost) implicit none real(dp), intent(in) :: chrom(:,:) integer , intent(in) :: i, npar !Local variables integer :: k real(dp) :: x(npar), cost do k = 1, npar x(k) = chrom(i,k) end do cost = x(1)**2 + x(2)**2 end function cost_func !***************************************************************************** !* Initial population !**************************************************************************** subroutine initial(chrom,npop,npar,pmx) implicit none real(dp) :: chrom(:,:),pmx(:,:),r integer :: i,j,npop,npar do i = 1, npop ! Assign a random number for each parameter. do j = 1, npar call random_number(r) chrom(i,j) = (pmx(2,j)-pmx(1,j))*r+pmx(1,j) end do ! Evaluate the cost function. chrom(i,npar+1) = cost_func(chrom,i,npar) end do end subroutine initial !***************************************************************************** !* Standard deviation for "nkeep" chromosomes. !**************************************************************************** function stndd(chrom,index,npar,nkeep,meanc) implicit none real(dp) :: chrom(:,:),std,meanc,stndd integer :: k,n,npar,nkeep,index(:) ; n=nkeep std=zero do k=1,n std = std + (chrom(index(k),npar+1)-meanc)**2 end do stndd = sqrt(std/real(n,dp)) end function stndd !**************************************************************************** !* heapsort: !* input : array ra(:) !* index indx(maxindx) !* integr desc give 1 if descending order is desired. !* for other values, ascending order will be assumed. !* output: array ra(:) !* index indx(maxindx) ordered as desired !**************************************************************************** subroutine heapsorti(ra,indx,maxindx,desc,npar) implicit none real(dp), intent(inout) :: ra(:,:) integer , intent(inout) :: indx(:) integer , intent(in) :: maxindx, desc, npar integer :: i,ir,j,m,nh,rra,nc nc = npar+1 if (maxindx < 2) return m=maxindx/2+1 ir=maxindx 10 continue if (m > 1) then m=m-1 rra=indx(m) else rra=indx(ir) indx(ir)=indx(1) ir=ir-1 if (ir == 1) then indx(1)=rra ! Reorder the array in descending order if (desc=1). if (desc == 1) then nh=maxindx/2 do i=1,nh rra=indx(i) indx(i)=indx(maxindx-(i-1)) indx(maxindx-(i-1))=rra end do endif return endif endif i=m j=m+m 20 if (j <= ir) then if (j < ir) then if (ra(indx(j),nc) < ra(indx(j+1),nc)) j=j+1 endif if (ra(rra,nc) < ra(indx(j),nc)) then indx(i) = indx(j) i = j j = j+j else j = ir+1 endif goto 20 endif indx(i)=rra goto 10 end subroutine heapsorti !**************************************************************************** !* Tournament selection !* 1. Select two chromosomes at random !* 2. Choose the one with smaller cost, return it as "tournament" !**************************************************************************** function tournament(chrom,index,nkeep,npar) implicit none real(dp) :: chrom(:,:), r integer :: nkeep, f1, f2, tournament, npar, index(:) call random_number(r) r = real(nkeep,dp)*r f1 = int(r)+1 call random_number(r) r = real(nkeep,dp)*r f2 = int(r)+1 if (chrom(index(f1),npar+1) < chrom(index(f2),npar+1)) then tournament = f1 else tournament = f2 endif end function tournament !**************************************************************************** !* Mating or crossover (producing an offspring) !* This produces an offspring "k" from the parents (np1,np2). !* blending for all the parameters. !**************************************************************************** subroutine crossover(chrom,index,k,np1,np2,npar) implicit none real(dp) :: chrom(:,:),beta integer :: i,k,np1,np2,npar,index(:) do i=1,npar call random_number(beta) chrom(index(k),i) = beta*chrom(index(np1),i) - (one-beta)*chrom(index(np2),i) end do do i=1,npar call random_number(beta) chrom(index(k+1),i) = - (one-beta)*chrom(index(np1),i) + beta*chrom(index(np2),i) end do end subroutine crossover !**************************************************************************** !* Mutation: (mrate = mutation rate) !**************************************************************************** subroutine mutation(chrom,index,npop,nkeep,npar,pmx) implicit none real(dp) :: chrom(:,:),pmx(:,:),r,mrate integer :: npop,nkeep,npar,i,imax,nm,nmp,index(:) mrate = 0.04_dp imax = int(mrate*npop*npar)+1 do i = 1, imax call random_number(r) r=real(npop,dp)*r nm=int(r)+1 if (nm==1) nm=nkeep+1 call random_number(r) r=real(npar,dp)*r if (r/=2.0d0) nmp=int(r)+1 call random_number(r) chrom(index(nm),nmp)=(pmx(2,nmp)-pmx(1,nmp))*r+pmx(1,nmp) end do end subroutine mutation !**************************************************************************** !***************************************************************************** end program genetic_algorithm