c	*************************************
c	* MULTIOBJECTIVE SIMULATED ANNEALING, MOSA-aJG*
c	*************************************
c	This is the present version of the code for 
c	Multi-Objective Simulated Annealing-aJG by SANKARARAO BODDUPALLI
c     The code for MOSA (developed by Sankararao B) is based on the work of 
c     Suppapitnarm (2000)
c     The ZDT4 problem is currently entered in the two subroutines, function1 and function2
c     Please keep the 'input.txt' file in the folder.
c     All rights reserved. This listing is for personal use.

	program mosa
	integer nt1,nt2,arnew,acc,arold,na,k1,k2,k3,nb,iter,b
	integer rstp,nh,ns,ns1,nacp
	real t1,t2,alow(20),ahigh(20),s(20),xn(20),xnew(20),un(20)
	real unew(20),x0(20),x,rb,ss,snew(20),c,ri,pjump
	real xnd(20000,20),xuncr(20000,20),xacc(60000,20),phi	!,xg(10000,20)
	integer tempe,jjump1,jjump2
	double precision fnd1(20000),ratio,sacc(60000,20)	!,fg1(10000),fg2(10000)
	double precision facc1(60000),facc2(60000),fold1,fold2
	double precision fnew1,fnew2,funcr1(20000),funcr2(20000)
	double precision fnd2(20000),fpara2(20000)
	double precision fdif1(20000),fdif2(20000),fpara1(20000) 	
      common/randvar/oldrand,jrand	
      real oldrand(55)
      integer jrand
	common/variables/nparam,randomseed 


c	open(unit=1,file='input.txt')

 	nparam = 10	
c	read(1,*)nparam				!nparam=no.of.parameters 	
	nobj = 2
	randomseed = 0.88876	
c	read(1,*)randomseed			!randomseed=value of randomseed 
	                            !between 0-1				
	iterations = 35000	
c	read(1,*)iterations			!iterations=no.of iterations 
	                            !to be performed      	
      ri = 0.9d0	
c	read(1,*)ri					!ri=value of the parameter ri
	rb = 0.9d0	
c	read(1,*)rb					!rb=value of the parameter rb
	b = 1	
c	read(1,*)b					!b=parameter for temperatuere(used 
	                            !in t=b*std 
	c = 0.4d0	
c	read(1,*)c					!c=parameter c(used in na=c*nt2)
	ss = 0.5d0	
c	read(1,*)ss					!ss=value of the step size		
	nt1 = 5000	
c	read(1,*)nt1				!nt1=no.of.iterations to be performed
	                            !before initial annealing schedule
	nt2 = 900	
c	read(1,*)nt2				!nt2=no=no.of.iterations to be 
	                            !performed for the next annealing schedule
	ns1 = 200	
c	read(1,*)ns1				!no of iterations to be performed before chang
								!ing stepsize
	pjump = 0.5					!jumping gene probability
c	close(1)



	write(*,*)'*****************input of SA****************'
	write(*,10)nparam
	write(*,11)randomseed
	write(*,12)iterations
	write(*,13)ri
	write(*,14)rb
	write(*,15)b
	write(*,16)c
	write(*,17)ss
	write(*,18)nt1
	write(*,19)nt2
	write(*,20)ns1
	write(*,21)pjump
10	format(2x,'The no of paramters are',i2)
11	format(2x,'The random seed is',f7.5)
12	format(2x,'The no of iterations are',i4)
13	format(2x,'The value of ri is',f7.5)
14	format(2x,'The value of rb is',f7.5)
15	format(2x,'The value of b is',i2)
16	format(2x,'The value of c is',f7.5)
17	format(2x,'The value of step size parameter is',f7.5) 
18	format(2x,'The value of nt1 is',i4)
19	format(2x,'The value of nt2 is',i4)
20	format(2x,'The value of ns1 is',i4)
21	format(2x,'The value of pjump is',f7.5)
	   		
      call randomize		

	arnew=0		 !arnew represents the no of uncrowded solutions taken
	             !after return to base 
	acc=1		 !acc represents no of accepted solutions
	arold=1		 !arold represents no of archiving solutions
	nacp=0
	ns=0
	phi=1.0	
	na=0		 !na represents no of acceptances after annealing scheduling
	nh=1		 !nh represents no of hundreds	
							
c	enter the bounds for deecision variables
	alow(1)=0.0
	ahigh(1)=1.0
      do j=2,nparam
	alow(j)=-5.0
	ahigh(j)=5.0
	end do

c	enter step size
	do j=1,nparam
      s(j)=0.5
	end do	
	
      do j=1,nparam
	sacc(1,j)=s(j)		  			
	end do

	do j=1,nparam
	call irnd(x,alow(j),ahigh(j))
	x0(j)=x
	end do

	do j=1,nparam	
	xn(j)=0.0
	xnew(j)=0.0
	un(j)=0.0
	unew(j)=0.0
	end do

c	enter the temperatures here
	t1=150.0 !t1 and t2 represents temperatures
	t2=150.0
	
	call function1(x0,funct1)
	fold1=funct1
	call function2(x0,funct2)
	fold2=funct2
	
c     fg1(1)=fold1		  !fg=the array contains all function values 
c	fg2(1)=fold2
	fnd1(1)=fold1		  !fnd=the array contains all non-dominated 
	                      !function values 
	fnd2(1)=fold2 
	facc1(1)=fold1		  !facc=the array contains all accepted 
	                      !function values 
	facc2(1)=fold2

	do j=1,nparam
c	xg(1,j)=x0(j)		  !xg=the array contains all values of x 
	xnd(1,j)=x0(j)		  !xnd=the array contains all non-dominating
	                      !values of x 
	xacc(1,j)=x0(j)		  !xacc=the array contains all accepted
	                      !values of x
	xn(j)=x0(j)
	end do
	
      !files for accepted,nondominated,uncrowded solutions
		
	open(unit= 3, file='nondom.txt')
      open(unit= 4, file='uncrowed.txt')	
	open(unit= 5, file='stepsize.txt')
	
	iter=1
c	initial perturbation to x............
222	call perturb(s,xn,xnew,un,unew,alow,ahigh)
	
      do i=1,nparam
      if((xnew(i).gt.ahigh(i)).or.(xnew(i).lt.alow(i)))then
	xnew(i) = alow(i) + (ahigh(i) - alow(i))*random()
	end if
	end do	

	If(iflip(pjump) .eq. 1) then
	jjump1 = nint(random()*(nparam-5))
	jjump2 = jjump1+5					!nint(random()*nparam)
	else
	jjump1 = jjump2
	goto  251
	endif
		 
c	if (jjump1 .gt. jjump2)	then 
c	tempe = jjump1
c	jjump1 = jjump2
c	jjump2 = tempe	 
c	endif
	do j1 = jjump1,jjump2
	xnew(j1) = alow(j1) + (ahigh(j1) - alow(j1))*random()
	end do     								  
                    
c	if((xnew(2).gt.ahigh(2)).or.(xnew(2).lt.alow(2)))
c	-	goto 222
c	if((xnew(1)**2+xnew(2)**2).gt.1.0) goto 222	  !constraint
251	iter=iter+1
	ns=ns+1
												   
	call function1(xnew,funct1)
	fnew1=funct1	
	call function2(xnew,funct2)
	fnew2=funct2	
			
c	fg1(iter)=fnew1
c	fg2(iter)=fnew2	
					
c	do j=1,nparam
c	xg(iter,j)=xnew(j)		 !all solutions
c	end do

	call archiving(fnew1,fnew2,fnd1,fnd2,xnd,k1,arold)			
	if(k1.eq.1)then		
		arold=arold+1		!no.of archive solutions
		acc=acc+1			!no.of accepted solutions
		nacp=nacp+1			!no.of accepted solutions for changing step size
		fnd1(arold)=fnew1
		fnd2(arold)=fnew2				
		facc1(acc)=fnew1
		facc2(acc)=fnew2
		do j=1,nparam
		xnd(arold,j)=xnew(j)	!archive solutions
		xacc(acc,j)=xnew(j)	!accepted solutions		
		xn(j)=xnew(j)
c		snew(j)=0.9*s(j)+ss*abs(unew(j)-un(j))		
c		if((snew(j).ge.1e-4).and.(snew(j).le.0.5)) s(j)=snew(j)						
		end do
		do j=1,nparam
		sacc(acc,j)=s(j)		  			
		end do
		if(iter.gt.nt1) na=na+1			
	else
	call probabilitytest(fnew1,fnew2,facc1,facc2,acc,t1,t2,k2)
		if(k2.eq.1)then			
			acc=acc+1
			nacp=nacp+1
			facc1(acc)=fnew1
			facc2(acc)=fnew2		
			do j=1,nparam
			xacc(acc,j)=xnew(j)			
			xn(j)=xnew(j)
c			snew(j)=0.9*s(j)+ss*abs(unew(j)-un(j))
c			if((snew(j).ge.1e-4).and.(snew(j).le.0.5)) s(j)=snew(j)
			end do
			do j=1,nparam
			sacc(acc,j)=s(j)		  			
			end do				
			if(iter.gt.nt1) na=na+1			
		end if
	end if
	
      if(ns.eq.ns1)then
	do i=1,nparam
	ratio=dfloat(nacp)/dfloat(ns)
      if(ratio.gt.0.6)then
      s(i)=s(i)*(1+2*(ratio-0.6)/0.4)
      else if(ratio.lt.0.4)then
      s(i)=s(i)/(1+2*((0.4-ratio)/0.4))
      end if
      if(s(i).gt.(ahigh(i)-alow(i)))then		
      s(i)=ahigh(i)-alow(i)	
      end if
	end do
	ns=0
	nacp=0
	write(5,204)iter
	write(5,*)s(1),s(2)
	end if
		  	
	if(iter.eq.nh*1000)then
	write(*,201)iter
	call writing(xnd,arold,fnd1,fnd2)
	write(*,*) acc

	open(unit= 2, file='accepted.txt')  	
	write(2,202)iter
	write(2,*)'Accepted x values are'
      do i=1,acc
 	write(2,*) xacc(i,1),xacc(i,2)
	end do
	write(2,*)'Accepted function values are'
	do i=1,acc
	write(2,*) facc1(i),facc2(i)
	end do
	close(2) 
			
	write(3,201)iter
	write(3,*)'Non-dominating x values are'
      do i=1,arold
 	write(3,*) xnd(i,1),xnd(i,2)
	end do
	write(3,*)'Non-dominating function values are'
	do i=1,arold
	write(3,*) fnd1(i),fnd2(i)
	end do
	
c	open(unit= 5, file='stepsize.txt')
c	write(5,204)iter                                
c     write(5,*)'The accepted values of stepsizes are'
c	do i=1,acc
c	write(5,*) sacc(acc,1),sacc(acc,2)
c	end do
c	close(5)
	
	if(nh.ge.6)then		
	write(4,203)iter
	write(4,*) 'Uncrowded x values are'
      do i=1,arnew
	write(4,*) xuncr(i,1),xuncr(i,2)
	end do
	write(4,*)'Uncrowded function values are'
	do i=1,arnew
	write(4,*) funcr1(i),funcr2(i)
	end do	
	end if
		
     	nh=nh+1	
c	pause
	end if

      if(iter.eq.nt1)then				
c		write(*,201)iter
201		format(2x,'The Archived solutions after',i5,'th iteration')	
202		format(2x,'The Accepted solutions after',i5,'th iteration')	
203		format(2x,'The Uncrowded solutions after',i5,'th iteration')
204		format(2x,'The accepted values of step sizes after',i5,'th
     +    iteration')

	    do i=1,arold
		fpara1(i)=0.0
		fpara2(i)=0.0
		end do
	
		do i=1,arold-1
		fdif1(i)=0.0
		fdif2(i)=0.0
		end do
     	
		do i=1,arold
		fpara1(i)=fnd1(i)
		end do

		call sort2(arold,fpara1)
		do i=1,arold-1
		fdif1=abs(fpara1(i+1)-fpara1(i))
		end do

		do i=1,arold
		fpara2(i)=fnd2(i)
		end do

		call sort1(arold,fpara2)
		do i=1,arold-1
		fdif2=abs(fpara2(i+1)-fpara2(i))
		end do


		call termination(k3,t1,t2,fdif1,fdif2,arold)
		if(k3.eq.1)then
		write(*,201)iter
		call writing(xnd,arold,fnd1,fnd2)
		goto 666
		end if

		call initannsche(t1,t2,facc1,facc2,acc,b)
		arnew=phi*arold
		call retbase(fnd1,fnd2,xnd,arold,arnew,xuncr,funcr1,funcr2)
	    rstp=nint(random()*(3-(0)))	 !arnew re-start point				
		do j=1,nparam
		xn(j)=xuncr(rstp,j)
		end do		

		write(4,203)iter
		write(4,*) 'Uncrowded x values are'
		do i=1,arnew
		write(4,*) xuncr(i,1),xuncr(i,2)
		end do
		write(4,*)'Uncrowded function values are'
		do i=1,arnew
		write(4,*) funcr1(i),funcr2(i)
		end do			 						
	end if
		
	if(iter.gt.nt1)then
c		no=no+1		!no=no.of.iterations to be performed for the next
		            !annealing schedule
		nb=nb+1		!nb=no.of.iterations to be performed for the next
		            !return to base	
c		if(nb.eq.int(rb*2*nt2))then
		if(nb.eq.200)then			
c			if(nb.lt.10) nb=10
c			if(nb.ge.10)then
c				arnew=int(ri*arold)
				arnew=arold
c				if(arnew.lt.6)then
c				arnew=6
c				write(*,201)iter
c				call writing(xnd,arold,fnd1,fnd2)
c				goto 666
c				end if				
				call retbase(fnd1,fnd2,xnd,arold,arnew,xuncr,
	&			funcr1,funcr2)
     				if(iter.le.55000)then												
				rstp=nint(random()*(3-(0)))	 !arnew re-start point								
				do j=1,nparam
				xn(j)=xuncr(rstp,j)
				end do
				else
				do j=1,nparam
				xn(j)=xuncr(3,j)
				end do
				end if
c				write(*,201)iter
c				call writing(xnd,arold,fnd1,fnd2)
c				pause
c				ri=0.9*ri
c				rb=0.9*rb
c				if(nb.eq.10) rb=0.9
				nb=0 		
c			else
c				write(*,201)iter
c				call writing(xnd,arold,fnd1,fnd2)
c				goto 666
c			end if
		end if

c		if(na.eq.int(c*nt2))then
		if(na.eq.nt2)then	
		
			do i=1,arold
			fpara1(i)=0.0
			fpara2(i)=0.0
			end do
		
			do i=1,arold-1
			fdif1(i)=0.0
			fdif2(i)=0.0
			end do

			do i=1,arold
			fpara1(i)=fnd1(i)
			end do

			call sort2(arold,fpara1)
			do i=1,arold-1
			fdif1(i)=abs(fpara1(i+1)-fpara1(i))
			end do

			do i=1,arold
			fpara2(i)=fnd2(i)
			end do

			call sort1(arold,fpara2)
			do i=1,arold-1
			fdif2(i)=abs(fpara2(i+1)-fpara2(i))
			end do
			
			call termination(k3,t1,t2,fdif1,fdif2,arold)
			if(k3.eq.1)then
			write(*,201)iter
			call writing(xnd,arold,fnd1,fnd2)
			goto 666
			end if

			call annsche(t1,t2,facc1,facc2,acc,na)

c			write(*,201)iter
c			call writing(xnd,arold,fnd1,fnd2)
c			pause
			na=0
c			no=0
		end if
				
c		if(no.eq.nt2)then	  		
	
c     			do i=1,arold
c			fpara1(i)=0.0
c			fpara2(i)=0.0
c			end do
			
c			do i=1,arold-1
c			fdif1(i)=0.0
c			fdif2(i)=0.0
c			end do
     		
c			do i=1,arold
c			fpara1(i)=fnd1(i)
c			end do

c			call sort2(arold,fpara1)
c			do i=1,arold-1
c			fdif1=abs(fpara1(i+1)-fpara1(i))
c			end do

c			do i=1,arold
c			fpara2(i)=fnd2(i)
c			end do

c			call sort1(arold,fpara2)
c			do i=1,arold-1
c			fdif2=abs(fpara2(i+1)-fpara2(i))
c			end do

c			call termination(k3,t1,t2,fdif1,fdif2,arold)
c			if(k3.eq.1)then
c			write(*,201)iter
c			call writing(xnd,arold,fnd1,fnd2)
c			goto 666
c			end if

c			call annsche(t1,t2,facc1,facc2,acc,na)

c			write(*,201)iter
c			call writing(xnd,arold,fnd1,fnd2)
c			pause
c			no=0
c			na=0
c		end if	
	end if						
      
	if(iter.lt.iterations)then
	goto 222
	else
	goto 666
c	write(*,201)iter
c	call writing(xnd,arold,fnd1,fnd2)
	end if		
	
666	close(3)
	close(4)
	close(5)
	
	stop
	end program mosa
	
      subroutine sort2(arold,fpara1)
      integer arold,i1,k1
      double precision temp,fpara1(arold)

      do k1=1,arold-1
        do i1=k1+1,arold
          if(fpara1(k1) .gt. fpara1(i1))then
             temp=fpara1(k1)             
             fpara1(k1)=fpara1(i1)             
             fpara1(i1)=temp             
          end if
        end do
      end do		  

	return
	end
	
      subroutine irnd(i,ilow,ihigh)	
	real ilow,ihigh,i
      if(ilow.ge.ihigh) then
         i = ilow
      else
         i = random()*(ihigh-ilow)+ilow
         if(i.gt.ihigh) i = ihigh
      endif      
      return
      end	
	
      subroutine retbase(fnd1,fnd2,xnd,arold,arnew,xuncr,funcr1,funcr2)
	common/variables/nparam,randomseed
	integer arold,arnew  
      double precision fnd1(arold),fnd2(arold),cub_len(arold)	
	real array2(20000,20),xnd(20000,20),xuncr(20000,20)	
	double precision fmax1,fmin1,fmax2,fmin2,array1(arold)
	double precision fpara(arold),length(arold),max,funcr2(arold)
	double precision array3(arold),array4(arold),funcr1(arold)	
	 	
      fmax1=maxval(fnd1)
	fmin1=minval(fnd1)
	fmax2=maxval(fnd2)
	fmin2=minval(fnd2)
	
	do i=1,arold      
	fnd1(i)=(fnd1(i)-fmin1)/(fmax1-fmin1)
	fnd2(i)=(fnd2(i)-fmin2)/(fmax2-fmin2)
	end do
	
      do i=1,arold
          fpara(i)=0.0		
          length(i)=0.0
		cub_len(i)=0.0		
      end do

	do i=1,arold
         fpara(i)=fnd1(i)
      end do

	do i=1,arold
		do j=1,nparam
		array2(i,j)=xnd(i,j)
		end do
	end do

	call sort(arold,fpara,array2)
	do i=1,arold
	array3(i)=fpara(i)	
	end do			 
      max = fpara(arold)      

	do i=1,arold
        if(i.eq.1 .or. i.eq.arold)then          
          length(i)=100*max
        else
          length(i)=abs(fpara(i+1)-fpara(i-1))
        end if
      end do

	do i=1,arold
       cub_len(i)=cub_len(i)+length(i)
      end do

	do i=1,arold
          fpara(i)=0.0          
          length(i)=0.0          
      end do

	do i=1,arold
         fpara(i)=fnd2(i)
      end do

	call sort1(arold,fpara)
	do i=1,arold
	array4(i)=fpara(i)
	end do			 
      max = fpara(arold)	

	do i=1,arold
        if(i.eq.1 .or. i.eq.arold)then          
          length(i)=100*max
        else
          length(i)=abs(fpara(i+1)-fpara(i-1))
        end if
      end do

	do i=1,arold
      cub_len(i)=cub_len(i)+length(i)
      end do

	do i =1,arold
	array1(i)=cub_len(i)
      end do	

	call sorting(array1,array2,array3,array4,arold) 

	do i=1,20000
	do j=1,nparam
	xuncr(i,j)=0.0
	end do
	end do

	if(arnew.le.arold)then
		do j=1,nparam
			xuncr(1,j)=array2(1,j)			
			do i=2,arnew-1
			xuncr(i,j)=array2(i,j)
			end do
			xuncr(arnew,j)=array2(arold,j)			
		end do
			funcr1(1)=array3(1)
			funcr2(1)=array4(1)
			do i=2,arnew-1			
			funcr1(i)=array3(i)
			funcr2(i)=array4(i)
			end do
			funcr1(arnew)=array3(arold)
			funcr2(arnew)=array4(arold)
			do i=1,arnew      
			funcr1(i)=fmin1+funcr1(i)*(fmax1-fmin1)
			funcr2(i)=fmin2+funcr2(i)*(fmax2-fmin2)
			end do
	else
		do i=1,arold
		funcr1(i)=array3(i)
		funcr2(i)=array4(i)		
		funcr1(i)=fmin1+funcr1(i)*(fmax1-fmin1)
		funcr2(i)=fmin2+funcr2(i)*(fmax2-fmin2)
		do j=1,nparam
		xuncr(i,j)=array2(i,j)
		end do
		end do
	end if

	do i=1,arold      
	fnd1(i)=fmin1+fnd1(i)*(fmax1-fmin1)
	fnd2(i)=fmin2+fnd2(i)*(fmax2-fmin2)
	end do
	
      do i=1,arold
		do j=1,nparam
		array2(i,j)=0.0
		end do
	end do
	
      do i =1,arold
	array1(i)=0.0
	array3(i)=0.0
	array4(i)=0.0
      end do
		
	return
	end

	subroutine sorting(array1,array2,array3,array4,arold)
	integer arold,i,j
      double precision array1(arold),array3(arold),array4(arold)
	double precision temp,temp1,temp2,temp3,temp4
	real array2(20000,20)

	do i=1,arold-1
       do j = i+1,arold
        if(array1(i) .le. array1(j)) then
         temp=array1(i)
	   temp1=array2(i,1)
	   temp2=array2(i,2)
	   temp3=array3(i)
	   temp4=array4(i)
	   array1(i)=array1(j)
	   array2(i,1)=array2(j,1)
	   array2(i,2)=array2(j,2)
	   array3(i)=array3(j)
	   array4(i)=array4(j)	   
         array1(j)=temp         
	   array2(j,1)=temp1
	   array2(j,2)=temp2         
	   array3(j)=temp3
	   array4(j)=temp4
        end if
       end do			 
      end do
	      
      return
      end

	subroutine sort(arold,fpara,array2)
      integer arold,i1,k1
      double precision temp,temp1,temp2,fpara(arold)
	real array2(20000,20)

      do k1=1,arold-1
        do i1=k1+1,arold
          if(fpara(k1) .gt. fpara(i1))then
             temp=fpara(k1)
		   temp1=array2(k1,1)
		   temp2=array2(k1,2)             
             fpara(k1)=fpara(i1)
		   array2(k1,1)=array2(i1,1)
		   array2(k1,2)=array2(i1,2)             
             fpara(i1)=temp             
		   array2(i1,1)=temp1
	       array2(i1,2)=temp2
          end if
        end do
      end do		  
	        
	return
      end

	subroutine sort1(arold,fpara)
      integer arold,i1,k1
      double precision temp,fpara(arold)

      do k1=1,arold-1
        do i1=k1+1,arold
          if(fpara(k1) .lt. fpara(i1))then
             temp=fpara(k1)             
             fpara(k1)=fpara(i1)             
             fpara(i1)=temp             
          end if
        end do
      end do		  
	        
		return
      end
	    
      subroutine initannsche(t1,t2,facc1,facc2,acc,b)
	real average1,average2,summation1,summation2,std1,std2	
	real sum1,sum2,t1,t2
	integer acc,b
	double precision facc1(acc),facc2(acc)

	sum1=0.0
	sum2=0.0
	
      do i=1,acc
	sum1=sum1+facc1(i)
	sum2=sum2+facc2(i)
	end do

	average1=sum1/acc
	average2=sum2/acc

	summation1=0.0
	summation2=0.0
	
     	do i=1,acc
	summation1=summation1+(facc1(i)-average1)**2
	summation2=summation2+(facc2(i)-average2)**2
	end do

	std1=(summation1/acc)**0.5
	std2=(summation2/acc)**0.5

	t1=b*std1
	t2=b*std2

	return
	end

	subroutine annsche(t1,t2,facc1,facc2,acc,na)
	real average1,average2,summation1,summation2,std1,std2	
	real sum1,sum2,t1,t2,alpha1,alpha2
	integer acc,na
	double precision facc1(acc),facc2(acc) 

	sum1=0.0
	sum2=0.0
	
      do i=acc-na+1,acc
	sum1=sum1+facc1(i)
	sum2=sum2+facc2(i)
	end do

	average1=sum1/na
	average2=sum2/na

	summation1=0.0
	summation2=0.0

	do i=acc-na+1,acc
	summation1=summation1+(facc1(i)-average1)**2
	summation2=summation2+(facc2(i)-average2)**2
	end do

	std1=(summation1/na)**0.5
	std2=(summation2/na)**0.5

	alpha1=max(0.5,exp((-0.7*t1)/std1))
	alpha2=max(0.5,exp((-0.7*t2)/std2))

	t1=alpha1*t1
	t2=alpha2*t2

	return
	end

	subroutine writing(xnd,j,fnd1,fnd2)
	common/variables/nparam,randomseed 
	integer j
	real xnd(20000,20)
	double precision fnd1(j),fnd2(j)	     	
      write(*,*) 'x values are'
      do i=1,j
	write(*,*) xnd(i,1),xnd(i,2)
	end do
	write(*,*)'function values are'
	do i=1,j
	write(*,*) fnd1(i),fnd2(i)
	end do
	
	return
	end
		
c     ***********put your 1st objectice function here **********
	subroutine function1(x,funct1)
	common/variables/nparam,randomseed
	real x(nparam),funct1
	if((x(1).le.0.01).and.(x(1).ge.0.0)) x(1)=0.0
	if((x(1).ge.0.99).and.(x(1).le.1.0)) x(1)=1.0
	do j=2,nparam
	if((x(j).le.0.01).and.(x(j).ge.0.0)) x(j)=0.0
	if((x(j).ge.-0.01).and.(x(j).le.0.0)) x(j)=0.0
	if((x(j).ge.4.99).and.(x(j).le.5.0)) x(j)=5.0
	if((x(j).le.-4.99).and.(x(j).ge.-5.0)) x(j)=-5.0
	end do				        
	funct1=x(1)		 
c	funct1=x(1)**2		 
	return
	end
c	 ***** put your 2nd objective function for optimization here *****
	subroutine function2(x,funct2)
	common/variables/nparam,randomseed
	real x(nparam),funct2,gx,sum       	
	
	if((x(1).le.0.01).and.(x(1).ge.0.0)) x(1)=0.0
	if((x(1).ge.0.99).and.(x(1).le.1.0)) x(1)=1.0
	do j=2,nparam
	if((x(j).le.0.01).and.(x(j).ge.0.0)) x(j)=0.0
	if((x(j).ge.-0.01).and.(x(j).le.0.0)) x(j)=0.0
	if((x(j).ge.4.99).and.(x(j).le.5.0)) x(j)=5.0
	if((x(j).le.-4.99).and.(x(j).ge.-5.0)) x(j)=-5.0
	end do
      sum=0.0
	do i=2,nparam
	sum=sum+(x(i)**2-10*cos(4*3.14159*x(i)))
	end do
	gx=1+10*(nparam-1)+sum	       	
	funct2=gx*(1-(x(1)/gx)**0.5)		 	
c	funct2=(x(1)-2)**2		 
	return
	end				 

c.....perturbation of xn.................
	subroutine perturb(s,xn,xnew,un,unew,alow,ahigh)
	common/variables/nparam,randomseed
	real s(nparam),xn(nparam),xnew(nparam),alow(nparam),ahigh(nparam)
	real un(nparam),unew(nparam)	
      do j=1,nparam
	xnew(j)=xn(j)+(-1+random()*(1-(-1)))*s(j)
c	un(j)=((xn(j)-alow(j))/(ahigh(j)-alow(j)))*2-1
c	unew(j)=un(j)+(-1+random()*(1-(-1)))*s(j)
c	xnew(j)=(((unew(j)+1)/2)*(ahigh(j)-alow(j)))+alow(j)
	end do						 
	return
	end    

      function modop(i,j)	
      modop = i - j*ifix(float(i)/float(j))
      return
      end	

	subroutine archiving(fnew1,fnew2,fnd1,fnd2,xnd,k1,arold)
	common/variables/nparam,randomseed
	integer k1,val,arold,count	!iter,
	real xnd(20000,20)
	double precision fnew1,fnew2,fnd1(arold),fnd2(arold)	
	k1=0
c	i=iter      
	do j=1,arold	 !	2
		val= value(fnd1(j),fnew1,fnd2(j),fnew2)
		if(val.eq.3)then
		k1=1
		else if(val.eq.1)then
		k1=1
		fnd1(j)=0.0
		fnd2(j)=0.0
		do k=1,nparam
		xnd(j,k)=0.0
		end do
		else
		k1=0
		goto 49
		end if
	end do
	
      count=0
      do i=1,arold 
		if((fnd1(i).ne.0.0).or.(fnd2(i).ne.0.0))then
		count=count+1
		fnd1(count)=fnd1(i)
     		fnd2(count)=fnd2(i)
		do k=1,nparam
		xnd(count,k)=xnd(i,k)
		end do
		end if
	end do

	arold=count
							 
49    return
      end	

	function value(fit11,fit12,fit21,fit22)
	double precision fit11,fit12,fit21,fit22	       
	if(((fit11 .ge. fit12) .and.((fit21 - fit22) .gt.
     1     1e-3)) .or.(((fit11 - fit12) .gt.1e-3) .and.
     1  (fit21.ge.fit22))) then
             value = 1
      else if(((fit11 .le. fit12).and.((fit22 - fit21).gt.
     1           1e-3)) .or.(((fit12 - fit11) .gt.1e-3) .and.
     1           (fit22 .ge. fit21))) then
             value= 2
      else 
             value=3	
      end if							 
      return
      end

	subroutine probabilitytest(f1,f2,fl1,fl2,j,t1,t2,k2)
	real p,p1,t1,t2
	integer k2,j
	double precision f1,f2,fl1(j),fl2(j)
	k2=0
c	if(((f1-fl1(j)).gt.(20.0*t1)).or.((f2-fl2(j)).gt.(20.0*t2)))then
c	k2=0
c	goto 453
c	end if	
c	if(((f1-fl1(j)).lt.(-20.0*t1)).or.((f2-fl2(j)).lt.(-20.0*t2)))then
c	k2=1
c	goto 453
c	end if
	if((t1.lt.4e-4).or.(t2.lt.4e-4))then
	k2=0
	goto 453
	end if
		
      p=(exp(-((f1-fl1(j))/t1))*exp(-((f2-fl2(j))/t2)))	
	p1 = random()												 
	if(p1.le.p) k2=1
												 
453	return
	end

	function iflip(prob)
      iflip = 0
      if ((prob.eq.1.0).or.(random().le.prob)) iflip = 1
      return
      end
     
	subroutine randomize
	common/variables/nparam,randomseed
      call warmup_random(randomseed)
      return
      end

	subroutine warmup_random(random_seed)
      real oldrand(55)
	common/randvar/oldrand,jrand            
      oldrand(55) = random_seed
      rand_new = 1.0e-9
      prev_rand = random_seed
      do 22 j1 = 1,54
         ii = modop(21*j1,55)
         oldrand(ii) = rand_new
         rand_new = prev_rand - rand_new
         if(rand_new.lt.0.0) rand_new = rand_new + 1.0
         prev_rand = oldrand(ii)
 22   continue
	call advance_random
      call advance_random
      call advance_random
      jrand = 0
      return
      end

	subroutine advance_random 
	real oldrand(55)
	common/randvar/oldrand,jrand      	
      do 13 j1 = 1,24
         rand_new = oldrand(j1) - oldrand(j1+31)
         if(rand_new.lt.0.0) rand_new = rand_new + 1.0
         oldrand(j1) = rand_new
 13   continue
      do 14 j1 = 25,55
         rand_new = oldrand(j1)-oldrand(j1-24)
         if(rand_new.lt.0.0)  rand_new = rand_new+1.0
         oldrand(j1) = rand_new
 14   continue
      return
      end

	function random()
      real oldrand(55)
	common/randvar/oldrand,jrand
      jrand = jrand + 1
      if(jrand.gt.55) then
         jrand = 1
         call advance_random
      endif	
      random = oldrand(jrand)     
	return
      end

	subroutine termination(k3,t1,t2,fdif1,fdif2,arold)
	common/variables/nparam,randomseed
	integer arold
	double precision fdif1(arold-1),fdif2(arold-1)
	double precision fmax1,fmax2
	integer k3
	real t1,t2
	fmax1=maxval(fdif1)	
	fmax2=maxval(fdif2)
	k3=0	
	if(((t1.lt.0.1e-5).or.(t2.lt.0.1e-5)).or.
     $((fmax1.le.1e-8).or.(fmax2.le.1e-8))) k3=1
	return
	end