Particle swarm optimization (입자 군집 최적화) [algorithm] by 바죠


입자 군집 최적화 [알고리듬]
Particle swarm  optimization (PSO)
http://en.wikipedia.org/wiki/Particle_swarm_optimization



http://en.wikipedia.org/wiki/Swarm_intelligence
__The_Swarm-original_paper.pdf
http://vimeo.com/17407010
http://www.mit.edu/~deweck/PDF_archive/3%20Refereed%20Conference/3_50_AIAA-2005-1897.pdf

 

입자 군집 최적화 방법은 컴퓨터 과학 분야에서 사용하는 일종의 광역 최적화 방법이다. 반복적 계산을 통하여 후보해들이 동시에 개선되게 하므로써 최종적으로 목적함수 최적화를 달성한다. 후보해들의 움직임은 특정한 수학 공식을 따르게 된다. 여러 개 후보해들을 동시에 취급하는 특징이 있다. 이점은 유전 알고리듬과 동일한 부분이다. 여러 후보해들 사이의 정보를 교환하는 특징이 있다. 하나의 후보해가 가지고 있는 정보와 무리 전체가 가지고 있는 정보늘 동시에 사용하는 특징이 있다. 유전 알고리듬에서 사용하는 변이와 교차라는 개념을 명시적으로 사용하지는 않는다. 대신에 기록된 자신의 역대 기록과 무리 내부에서 공개된 역대 기록을 동시에 참고하고 직접 사용한다.  



아주 간단한 광역 최적화(global optimization) 알고리듬이다.  아마도 가장 간단한 광역 최적화 알고리듬이다. 유전 알고리듬(genetic  algorithm)과 비교할 때 매우 단순한 알고리듬이라고 볼 수 있다. 유전 알고리듬을 이해하고 있다면 금방 PSO 를 이해할 수 있다.  그 반대도 마찬가지이다. 하지만, 언급한 두 가지 알고리듬들 사이에는 개념적으로, 근본적으로 많은 차이가 있다. 실제 연습에서 유전 알고리듬(교차, 변이) 보다 먼저 입자 군집 최적화 알고리듬을 구현해 보는 것이 결코 나쁘지 않다. 왜나 하면, 입자 군집 최적화 알고리듬이 보다 더 간단하기 때문이다. 새, 물고기와 같이 집단행동을 하는 생물들의 집단행동 특성에서 광역 최적화의 아이디어를 얻었다. PSO와 유전 알고리듬의 공통점은 두 방법 모두 한 무리의 해들을 가정하는 것에서 찾을 수 있다. 두 방법 모두 병렬화가 용이하다. 병렬 효율성도 매우 뛰어나다. 함수의 일차미분값들을 필요로 하지 않는다. 따라서, 일차 미분값들이 알려지지 않은 매우 일반적인 복잡한 함수 최적화에 쉽게 적용될 수 있다.



https://en.wikipedia.org/wiki/James_Kennedy_(social_psychologist)
https://en.wikipedia.org/wiki/Russell_C._Eberhart


가장 중요한 점은 입자(새 한 마리, paticle, 후보해)의 운동은 무리, 군집(swarm)의 정보를 적극적으로 이용한다는 점이다. 사실상 정보교환을 통해서 얻어진 군집(무리)의 결정을  매우 중요하게 고려하는 것이다. 물론, 새 한 마리도 무리의 결정을 위해서 정보를 공유한다. 새 한 마리 입장에서, 자신의 역대  기록을 제공한다. 또한, 군집 전체의 역대 정보를 활용한다. 여기서 역대 기록이라고 함은 입자 하나가 시간이 지나면서 얻어낸 최소 목적  함수값을 제공하는 위치를 지칭하는 것이다. 군집 전체에 동일한 조건을 요구할 수 있다. 군집 전체가 지금까지 얻어낸 최소 목적 함수 값을 주는  입자의 위치를 보관할 필요가 있다. 유전 알고리듬에서와 마찬가지로 다수의 해(입자의 위치, 새 한마리의 위치)들을 동시에 취급한다.  또한, 확률론적 접근법이다. 결국, 많은 계산을 수반할 수밖에 없다.  얻어진 해가 절대적인 광역 최적화의 끝이라는 것을 알 수 없다. 


입자의 위치가 해가 된다. 따라서 자연스럽게 다차원 벡터량으로 해를 표현할 수 있다.  비트열(bit string)같은 것을 사용하지 않는다. 실수형 숫자로 표현된 벡터량 표현으로 입자의 위치벡터(후보해)를 표현한다.

x ---->  x' 처럼 시도 해를 만들어 낼 때 랜덤 넘버(막수, 무작위 숫자, 0에서 1사이의 임의의 숫자)도 활용하여 계산한다. 목적 함수 f(x)를  최소화한다.  여기에서는  목적 함수 최소화나 최대화나 같은 의미로 받아 드릴 수 있다. -1 이라는 숫자를 목적 함수에 곱하고,  새롭게 목적 함수를 정의하면 결국  목적 함수 최대화와 최소화는 실제 계산에서는 사실상 같은 개념의 것이 된다. x는 해를 나타내는 벡터량이다. 다수의 해들({x})을 동시에 취급한다.  각각의 해 x 는 기본적으로 변환된다. x + v ---> x 형식으로  변환되어진다. 여기에서  v 는 정보공유로 부터 나오는 계산 결과물이다. 통상 입자의 속도라고 부른다.무리로 부터의 정보를 제공받게 된다.


그림 참조:

https://www.google.co.kr/search?q=%EA%B0%80%EC%B0%BD%EC%98%A4%EB%A6%AC+%EA%B5%B0%EB%AC%B4&newwindow=1&biw=1462&bih=2407&tbm=isch&imgil=ZcXo9G-1BhLllM%253A%253BKoCKHGOdXA6eDM%253Bhttp%25253A%25252F%25252Fwww.namdokorea.com%25252Fen%25252Ftourism%25252F01%25252Fview02.jsp%25253Ftype%2525253D%25252526kind%2525253D16%25252526tour_id%2525253D335%25252526page%2525253D&source=iu&pf=m&fir=ZcXo9G-1BhLllM%253A%252CKoCKHGOdXA6eDM%252C_&usg=__zjWKzJ8AJVX-9YUqgDQVmzTqCO8%3D&ved=0ahUKEwiqx9Gur8fMAhXKq5QKHagVCHQQyjcIKw&ei=sY8tV-qqF8rX0gSoq6CgBw#imgrc=tKlqZq19YdITYM%3A

사회적 무리의 운동에서 아이디어를 얻고 광역 최적화 방편으로 개발된 알고리듬이다.  새떼, 고기떼 집단행동 특성에서 최적화 문제 풀이 알고리듬의 아이디어를 취한 것이라고 볼 수 있다.  각 개체의 역사와 무리의 역사를 모니터링 하게 되어 있다. 개인 역대 최적점, 무리 역대 최적점같은 기록들을 소중하게 생각한다. 이들을 지속적으로  컴퓨터 메모리에 저장한다. 개체의 특성이 유전된다는 개념을 활용하지 않는다.  유전 알고리듬(GA)의 용어가 전혀 사용되지 않는다.


 

 



사회적 무리의 운동에서 아이디어를 얻고 광역 최적화 방편으로 개발된 알고리듬이다.  새떼, 고기떼 집단행동 특성에서 최적화 문제 풀이 알고리듬의 아이디어를 취한 것이라고 볼 수 있다.  각 개체의 역사와 무리의 역사를 모니터링 하게 되어 있다. 개인 역대 최적점, 무리 역대 최적점같은 기록들을 소중하게 생각한다. 이들을 지속적으로  컴퓨터 메모리에 저장한다. 개체의 특성이 유전된다는 개념을 활용하지 않는다.  유전 알고리듬(GA)의 용어가 전혀 사용되지 않는다.

 

PSO는 유전 알고리듬과 전혀 다른 개념에 기반 한다. 개체들의 운동 특성과 개체간의 정보교환에 집중한다.  유전 알고리듬 (genetic algorithm; GA)에서 중요한  개념인 crossover, mutation이라는개념이 활용되지 않는다는 점에 주의해야만 한다.  또한, 선택, 대치의 개념도 없다. 우월한 유전자를 찾아내고 장려하는 선택이 없다. 또한 입자는 결코 다른 입자들에 의해서 대체되지 않는다. 다만, 변화 할뿐이다. 표면적으로는 그렇다. 하지만, 나중에 실제 계산 작업 단계로 가면 유사성을 얻어 낼 수도 있다.

비록 수치적으로 같은 개념일 수  있지만, 실제로 도입되는 단계에서는 결코, crossover, mutation이라는 개념을  동원하지 않는다.  대신에, 무리들의 운동에서 개체의 운동 특성을 존중한다.  또한 무리들 내부의 공통의 기록과 정보를 반드시 공유하게 된다.  개인과 집단  전체 사이의 관계라는 관점에서 이해해야 하는 PSO 알고리듬이다.

집단의 수행 특성을 기준으로 새로운 해들을 찾아 나가게 된다.  개인의 기록, 최적 기록과 집단의  기록, 최적 기록을각각 기록해 보관하는 특징이 있다.  PSO 에서는 개체간의 정보 공유가 키워드가 될 수 있다.  유전체와 진화를 다루는 것은 유전 알고리듬이다. 키워드가 사뭇 다르다.  결국, 유전 알고리듬의 핵심 사항이고 문제마다 새롭게 정의해야 하는 crossover,  mutation을 추가로 다룰 필요가 없어진다. selection pressure 도 다루지 않아도 됩니다.  아무튼, 무지하게 간단해집니다.   여러 해를 동시에 취급하는 방법들 중에서 아마도 가장 간단한 알고리듬이 아닐까? 다시  말해서, 유전적 특징에 기인한 이론이 아니다. 오히려,  행동의 특성과 무리 안의 정보 공유 특성에 기반한  알고리듬이다. 유전 알고리듬은 염색체 정보가 특성을 결정할 것이라고 믿는 반면,  PSO에서는 보여준 개체들 특성에 보다 더  집중한다. 개체와 군집(무리)의 특성들 사이에서 새로운 특성을 찾아내게 된다.  개체와 무리의 상호작용을 직접적으로 다룬다. 군집(무리) 사회에서의 정보 공유가 기본적인 가정으로 여겨진다.  즉, 유전 알고리듬에서 활용하는 개념인 선택이 여기에서도 존재한다.  사실상, 통신을 통해서,  어떠한 개체가 현재까지 이룩한 최고의 상황(목적 함수 최소값을 주는 개체의 위치)을 모든 개체가 다 공유하게 되어 있다.

개체 최적 상황(목적 함수 값이 가장 작은 경우에 해당하는 위치)도 각 개체는 기억을 하고 있다. 여전히 집단 ---> 집단의 방식으로 새로운 개체(개체의 위치)를 찾는다.

{x}  --->  {x'}결국, 새로운 개체(개체의 위치)를 만들어내는 방법이 다른 것이다.  새로운 개체(개체의 위치)를 동시에 취급하고 다루는 방법에서, 알고리듬 기본 이론상으로 가장 밑바탕이 되는 개념에는, 유전 알고리듬과 큰 차이가 없다고 볼 수도 있다. 다만, PSO에서는 mutation, crossover를 활용하지 않는다는 명백한 차이점은 존재한다.  입자 위치: 후보 해, candidate solution 여러 개의 particle들을 동시에 다룬다. 각자  입자의 위치를 가지고 있다. 이 위치가 계속해서 바뀔 것이다. 한 입자에 국한해서 찾은 최고의solution, 또한 다른 입자들까지 포함해서 찾은 최고의 solution등을 계속 모니터링 할 필요가 있다. 이 알고리듬에서 매우 중요한 것은 아래의 두 가지이다.계속해서 반복적으로  새로운 입자의 위치를 찾는다.  그 와중에 아래의 두 가지 양들을 계속해서 찾아주어야 한다. particle's best known position, swarm's best  known position

군집, 무리 역대 최적 위치(swarm's best  known position) 입자(particle) 위치가 있고  이들의 집단인  군집(swarm)이 있다. 역대 입자의 최적 위치를 정의할 수 있다.  Xbest, 모든 입자들로부터 Xbest를 제출받아서 최고중의  최고 위치를 정의할 수 있다. Gbest라고 부른다. 역대 군집 최적 위치라고 부를 수 있다.

이 두 가지 벡터량을 항상  정의하고 보유 할 수 있어야 한다. 이들로 부터 새로운 입자 위치가 찾아지게 된다. 이들을 참조하는 것이 유전 알고리듬에서 이야기하는 mutation과 유사한 기능을 수행하게 된다.  simulated  annealing(SA), parallel tempering(PT) 방법과 마찬가지로 기본적으로 엄청난 함수 계산을수행할 수 있어야만 적용 가능한 알고리듬이다.

 

분명히 PSO가  SA, PT, GA 보다 더 간단한 알고리듬처럼 보인다. 부가적으로 만들어야 하는 컴퓨터 계산 루틴들의 숫자가 줄어든다. 예를 들어, 유전  알고리듬에서 필요한 mutation, crossover, selection pressure, selection, replacement, generation, elitism  등과 연관된 부가적인 루틴들이 필요 없다.  만약 유전 알고리듬과 직접 비교한다면 많은 수의 루틴들이 필요 없다.  population size (군집의 크기, 군집을 이루는 개체의 수) 는 대략 적으로10+[sqrt(2 D)]정도로 잡는다.  여기서 D는 결정되어야 할 숫자들의 수를 나타낸다.

 

실제 응용에서는 국소 최소화 알고리듬을 활용한다. BFGS, Nelder-Mead 알고리듬을 활용한다. 전자는 목적함수의 해석적 미분이 있을 때 사용하는 국소 알고리듬이다. 후자는 목적함수의 해석적 미분이 없을 때 사용하는 알고리듬이다.  다시 말해서, 국소 최소화된 해만을 취급하는 PSO 방법이 더 효율적이다. 이러한 경향은 유전 알고리듬에서도 확인된다.


PSO 에서는 최고 성능 해의 영향력이 절대적이다. 늘 영향력을 행사하고 있다. 이 부분이 GA(유전 알고리듬)과 차이점이다. 따라서, PSO 에서는 현재 최적화되어 있는 해를 국소 최적화시키는 작업을 추가적으로 진행할 수 있다.  이 부분은 광역최적화 작업에서 중요한 차이점을 유발할 수 있다. PSO 또는 GA 모두 국소 최적화 작업을 수행해야만 효율적으로 광역최적화를 이룰 수 있다.

 

유전 알고리듬(GA)과 PSO는 모두 다 후보해들의 차이에 대한 언급이 전혀 없다.

다시 말해서, 고려하고 있는 유한한 갯수의 후보해들에서 지속적으로 다양성을 확보할 수 있는 방법이 없다.

이러한  문제점은 GA/PSO 기반 광역 최적화의 어려움을 낳게 된다.

이러한 약점을 해결할 수 있는 방법이 CSA 방법이다.

http://incredible.egloos.com/6208916

http://webzine.kps.or.kr/contents/data/webzine/webzine/15211841581.pdf


 

 병렬 조질, 풀리 시늉, 유전 알고리듬:
http://incredible.egloos.com/4784590
http://incredible.egloos.com/4784592
http://incredible.egloos.com/4856511



--------------------------------------------


!234567890
!      Written by In-Ho Lee, KRISS, July 18, 2018.
       module pso
       implicit none
       private
       save
       integer ndf,npop,nlocal,mpsoupdate
       integer*8 nfc
       real*8 xn1,xn2,xk
       real*8, allocatable :: xxx(:,:),yyy(:,:),vvv(:,:),arrin(:),energy(:)
       real*8, allocatable :: xlw(:),xup(:),x(:),y(:),epbest(:),pbest(:,:),gbest(:)
       integer, allocatable :: indx(:)

       public :: pso_initial,pso_final,update

       contains
!234567890
!      Written by In-Ho Lee, KRISS, July 18, 2018.
       subroutine pso_initial(ndf0,npop0,mpsoupdate0)
       implicit none
       integer ndf0,npop0,mpsoupdate0
       real*8 tmp,test0
       integer i,j,j1

       ndf=ndf0
       npop=npop0
       mpsoupdate=mpsoupdate0
       nlocal=10000
       allocate(x(ndf),y(ndf))
       allocate(xlw(ndf),xup(ndf))
       allocate(xxx(ndf,npop))
       allocate(yyy(ndf,npop))
       allocate(vvv(ndf,npop))
       allocate(gbest(ndf))
       allocate(pbest(ndf,npop))
       allocate(epbest(npop))
       allocate(energy(npop))
       allocate(arrin(npop))
       allocate(indx(npop))
       do i=1,ndf
       xlw(i)=-3.d0
       xup(i)= 3.d0
       enddo
       call random_seed()
       do j=1,npop
       do i=1,ndf
       call random_number(tmp)
       xxx(i,j)=xlw(i)+tmp*(xup(i)-xlw(i))
       call random_number(tmp)
       vvv(i,j)=-(xup(i)-xlw(i))+tmp*((xup(i)-xlw(i))+(xup(i)-xlw(i)))
       enddo
       enddo
       nfc=0
       do j=1,npop
       epbest(j)=objft(xxx(1,j))
       pbest(:,j)=xxx(:,j)
       enddo
       test0=minval(epbest)
       do j=1,npop
       if(abs(test0-epbest(j)) < 1.d-9)then
       gbest(:)=pbest(:,j)
                                       endif
       enddo
       arrin=energy
       call sortnr(npop,arrin,indx)
       do j=1,npop
       j1=indx(j) ; yyy(:,j)=xxx(:,j1) ; energy(j)=arrin(j1)
       enddo
       xn1=2.05d0 ; xn2=2.05d0
       tmp=xn1+xn2
       xk=2.d0/abs(2.d0-tmp-sqrt(tmp**2-4.d0*tmp))
       write(6,'(3e18.8,1x,a10)') xn1,xn2,xk,'xn1,xn2,xk'
       end subroutine pso_initial
!234567890
!      Written by In-Ho Lee, KRISS, July 18, 2018.
       subroutine pso_final()
       implicit none

       deallocate(x,y)
       deallocate(xlw,xup)
       deallocate(xxx,yyy,vvv)
       deallocate(epbest,pbest,gbest)
       deallocate(energy)
       deallocate(arrin)
       deallocate(indx)
       end subroutine pso_final
!234567890
!      Written by In-Ho Lee, KRISS, July 18, 2018.
       subroutine update()
       implicit none
       integer k,j,j1

       do k=1,mpsoupdate
       call swarmmove()

       arrin=epbest
       call sortnr(npop,arrin,indx)
       do j=1,npop
       j1=indx(j) ; energy(j)=arrin(j1) ; yyy(:,j)=pbest(:,j1)
       enddo

       call onedprint10(energy,npop)
       call onedprint10(yyy(1,1),ndf)
       enddo
       end subroutine update
!234567890
!      Written by In-Ho Lee, KRISS, July 18, 2018.
       subroutine swarmmove()
       implicit none
       integer i,j,klocal
       real*8 test0,tmp,tmq,tmr,zzr,zzi
       logical loutside

       do klocal=1,nlocal
       do j=1,npop
       do i=1,ndf
       call random_number(zzr) ; call random_number(zzi)
       vvv(i,j)=xk*vvv(i,j)+xn1*zzr*(pbest(i,j)-xxx(i,j))+xn2*zzi*(gbest(i)-xxx(i,j))
       xxx(i,j)=xxx(i,j)+vvv(i,j)
       loutside=.false.
       if(xxx(i,j) > xup(i) .or. xxx(i,j) < xlw(i)) loutside=.true.
       if(loutside)then
       call random_number(tmp) ; xxx(i,j)=xlw(i)+tmp*(xup(i)-xlw(i))
                   endif
       enddo
       test0=objft(xxx(1,j))
       if(test0 < epbest(j))then
       epbest(j)=test0
       pbest(:,j)=xxx(:,j)
                            endif
       enddo
       test0=minval(epbest)
       do j=1,npop
       if(abs(test0-epbest(j)) < 1.d-9)then
       gbest(:)=pbest(:,j)
                                       endif
       enddo
       enddo
       end subroutine swarmmove
!234567890
!      Written by In-Ho Lee, KRISS, July 18, 2018.
       real*8 function objft(x)
       implicit none
       real*8 x(ndf)
       real*8 pi
       integer i

       objft=0.d0
       pi=4.d0*atan(1.d0)
       objft=10.d0*ndf
       do i=1,ndf
       objft=objft+x(i)**2-10.d0*cos(2.d0*pi*x(i))
       enddo
       nfc=nfc+1
       end function objft
       end module pso
!234567890
!      Written by In-Ho Lee, KRISS, July 18, 2018.
       program pso_test
       USE pso, ONLY : pso_initial,pso_final,update
       implicit none
       integer ndf0,npop0,mpsoupdate0

!      call timestamp()
       ndf0=10
       npop0=50
       mpsoupdate0=10000
       call pso_initial(ndf0,npop0,mpsoupdate0)
       call update()
       call pso_final()
!      call timestamp()
       end program pso_test
!234567890



--------------------------------------------

!      Rastrigin function
       if(kase == 1)then
       pi=4.d0*atan(1.d0)
       fn=10.d0*nvar
       do i=1,nvar
       fn=fn+x(i)**2-10.d0*cos(2.d0*pi*x(i))
       enddo
                    endif
!      Ackley function
       if(kase == 2)then
       pi=4.d0*atan(1.d0)
       aa=20.d0 ; bb=0.2d0 ; cc=2.d0*pi ; tmq=0.d0 ; tmr=0.d0
       do i=1,nvar
       tmq=tmq+x(i)**2 ; tmr=tmr+cos(cc*x(i))
       enddo
       tmq=tmq/nvar
       tmq=-bb*sqrt(tmq)
       tmr=tmr/nvar
       if(tmq < -80.d0) tmq=-80.d0 ; if(tmq > 80.d0) tmq=80.d0
       if(tmr < -80.d0) tmr=-80.d0 ; if(tmr > 80.d0) tmr=80.d0
       fn=exp(1.d0)+aa-aa*exp(tmq)-exp(tmr)
                    endif
!      Eggholder function
       if(kase == 3)then
       fn=0.d0
       do i=1,nvar-1
       fn=fn-(x(i+1)+47.d0)*sin(sqrt(abs(0.5d0*x(i)+x(i+1)+47.d0))) &
                  -x(i)*sin(sqrt(abs(x(i)-x(i+1)-47.d0)))
       enddo
                    endif
!      Rosenbrock function
       if(kase == 4)then
       fn=0.d0
       do i=1,nvar-1
       fn=fn+(100.d0*(x(i+1)-x(i)**2)**2+(1.d0-x(i))**2)
       enddo
                    endif
!      Sphere function
       if(kase == 5)then
       fn=0.d0
       do i=1,nvar
       fn=fn+x(i)**2
       enddo
                    endif
!      Styblinski-Tang function
       if(kase == 6)then
       fn=0.d0
       do i=1,nvar
       fn=fn+x(i)**4-16.d0*x(i)**2+5.d0*x(i)
       enddo
       fn=fn/2.d0
                    endif
!      Griewank function
       if(kase == 7)then
       tmr=1.d0
       do i=1,nvar
       tmr=tmr*cos(x(i)/sqrt(dble(i)))
       enddo
       fn=-tmr+1.d0
       do i=1,nvar
       fn=fn+x(i)**2/4000.d0
       enddo
                    endif
!      Zakharov function
       if(kase == 8)then
       fn=0.d0
       tmr=0.d0
       do i=1,nvar
       fn=fn+x(i)**2
       tmr=tmr+0.5d0*dble(i)*x(i)
       enddo
       fn=fn+tmr**2+tmr**4
                    endif
!      Schwefel function
       if(kase == 9)then
       fn=nvar*418.9829d0
       do i=1,nvar
       fn=fn-x(i)*sin(abs(x(i)))
       enddo
                    endif
!      Levy function
       if(kase == 10)then
       pi=4.d0*atan(1.d0)
       fn=(sin(pi*(1.d0+(x(i)-1.d0)/4.d0)))**2
       do i=1,nvar
       tmq=1.d0+(x(i)-1.d0)/4.d0
       fn=fn+(tmq-1.d0)**2 *(1.d0+10.d0*(sin(pi*tmq+1.d0))**2)
       enddo
       tmq=1.d0+(x(nvar)-1.d0)/4.d0
       fn=fn+(tmq-1.d0)**2 *(1.d0+(sin(2.d0*pi*tmq))**2)
                     endif
!      Powell function
       if(kase == 11)then
       fn=0.d0
       do i=1,nvar/4
       tmq=x(4*i-3)+10.d0*x(4*i-2)
       fn=fn+tmq**2
       tmq=x(4*i-1)-x(4*i)
       fn=fn+5.d0*tmq**2
       tmq=x(4*i-2)-2.d0*x(4*i-1)
       fn=fn+tmq**4
       tmq=x(4*i-3)-x(4*i)
       fn=fn+10.d0*tmq**4
       enddo
                     endif



       if(kase == 1)then
       do i=1,nvar
       xl(i)=-5.12d0-1.d-4
       xu(i)= 5.12d0+1.d-4
       enddo
                    endif
       if(kase == 2)then
       do i=1,nvar
       xl(i)=-5.d0-1.d-4
       xu(i)= 5.d0+1.d-4
       enddo
                    endif
       if(kase == 3)then
       do i=1,nvar
       xl(i)=-512.d0-1.d-4
       xu(i)= 512.d0+1.d-4
       enddo
                    endif
       if(kase == 4)then
       do i=1,nvar
       xl(i)=-1.d10-1.d-4
       xu(i)= 1.d10+1.d-4
       enddo
                    endif
       if(kase == 5)then
       do i=1,nvar
       xl(i)=-1.d10-1.d-4
       xu(i)= 1.d10+1.d-4
       enddo
                    endif
       if(kase == 6)then
       do i=1,nvar
       xl(i)=-5.d0-1.d-4
       xu(i)= 5.d0+1.d-4
       enddo
                    endif
       if(kase == 7)then
       do i=1,nvar
       xl(i)=-600.d0-1.d-4
       xu(i)= 600.d0+1.d-4
       enddo
                    endif
       if(kase == 8)then
       do i=1,nvar
       xl(i)=-5.d0-1.d-4
       xu(i)=10.d0+1.d-4
       enddo
                    endif
       if(kase == 9)then
       do i=1,nvar
       xl(i)=-500.d0-1.d-4
       xu(i)=500.d0+1.d-4
       enddo
                    endif
       if(kase == 10)then
       do i=1,nvar
       xl(i)=-10.d0-1.d-4
       xu(i)=10.d0+1.d-4
       enddo
                     endif
       if(kase == 11)then
       do i=1,nvar
       xl(i)=-4.d0-1.d-4
       xu(i)=5.d0+1.d-4
       enddo
                     endif

--------------------------------------------


       module modpostest
       implicit none
       private
       save
       integer nvar00
       real*8 amp
       real*8, allocatable :: qosi0(:),posi(:,:)
       public :: qosi0,posi,amp,nvar00,fn

       contains
       real*8 function fn(x)
       implicit none
       real*8 x(nvar00)
       integer i

       fn=0.d0

       do i=1,nvar00
       fn=fn+(dble(i)-x(i))**2
       enddo
       return
       end function fn
       end module modpostest
!234567890
!      Written by In-Ho Lee, KRISS, August 1, 2018.
       module pso
       implicit none
       private
       save
       integer ndf,npop,mpsoupdate,nfactor,nsteps
       real*8 xn1,xn2,xk,xlocal
       real*4, allocatable :: vvv(:,:)
       real*8, allocatable :: epbest(:),pbest(:,:),gbest(:)
       real*8, allocatable :: omega(:),cc11(:),cc22(:)
       public :: pso_initial,pso_final,pso_do
       contains
!234567890
!      Written by In-Ho Lee, KRISS, August 1, 2018.
       subroutine pso_initial(ndf0,npop0,mpsoupdate0,nfactor0)
       USE modpostest, ONLY : nvar00,amp,qosi0,posi,fn
       implicit none
       integer ndf0,npop0,mpsoupdate0,nfactor0
       real*8 tmp,test0
       integer i,j

       nvar00=ndf0 ; amp=2.d0*nvar00
       allocate(qosi0(nvar00)) ; allocate(posi(nvar00,npop0))
       call random_number(qosi0) ; call random_number(posi)

       ndf=ndf0 ; npop=npop0 ; mpsoupdate=mpsoupdate0 ; nfactor=nfactor0
       nsteps=1000
       allocate(vvv(ndf,npop),pbest(ndf,npop),gbest(ndf),epbest(npop))

       xlocal=0.1d0
       call random_number(tmp)
       call random_number(vvv) ; vvv=(vvv-0.5d0)*2.d0 ; vvv=vvv*(xlocal)*tmp
       do j=1,npop
       posi(:,j)=qosi0(:)
       if(j /= 1) posi(:,j)=posi(:,j)+vvv(:,j)
       epbest(j)=fn(posi(1,j)) ; pbest(:,j)=posi(:,j)
       enddo
       test0=minval(epbest)
       do j=1,npop
       if(abs(test0-epbest(j)) < 1.d-9)then
       gbest(:)=pbest(:,j)
                                       endif
       enddo
       xn1=2.05d0 ; xn2=2.05d0
       tmp=xn1+xn2 ; xk=2.d0/abs(2.d0-tmp-sqrt(tmp**2-4.d0*tmp))
       allocate(omega(npop),cc11(npop),cc22(npop))
       do j=1,npop
       call random_number(tmp) ; omega(j)=0.4d0+(1.4d0-0.4d0)*tmp
       call random_number(tmp) ; cc11(j)=1.5d0+(2.0d0-1.5d0)*tmp
       call random_number(tmp) ; cc22(j)=2.0d0+(2.5d0-2.0d0)*tmp
       enddo
!      omega=xk ; cc11=xn1 ; cc22=xn2
!      omega=0.5d0 ; cc11=1.5d0 ; cc22=1.5d0
!      omega=1.0d0 ; cc11=2.0d0 ; cc22=2.0d0
       end subroutine pso_initial
!234567890
!      Written by In-Ho Lee, KRISS, August 1, 2018.
       subroutine pso_final(test)
       USE modpostest, ONLY : qosi0,posi
       implicit none
       real*8 test

       test=minval(epbest) ; qosi0=gbest
       write(6,*) test
       write(6,*) qosi0
       deallocate(qosi0) ; deallocate(posi)
       if(allocated(vvv)) deallocate(vvv)
       if(allocated(epbest)) deallocate(epbest)
       if(allocated(pbest)) deallocate(pbest)
       if(allocated(gbest)) deallocate(gbest)
       if(allocated(omega)) deallocate(omega)
       if(allocated(cc11)) deallocate(cc11)
       if(allocated(cc22)) deallocate(cc22)
       end subroutine pso_final
!234567890
!      Written by In-Ho Lee, KRISS, August 1, 2018.
       subroutine pso_do()
       USE modpostest, ONLY : amp,posi
       implicit none
       integer k,j,i
       real*8 tmp

       call random_number(vvv) ; vvv=(vvv-0.5d0)*2.d0*(xlocal) ; vvv=vvv*amp
       do k=1,mpsoupdate
       call swarmmove()
       if(mpsoupdate /= 1)then
       do j=1,npop
       do i=1,ndf
       call random_number(tmp) ; posi(i,j)=posi(i,j)+2.d0*(tmp-0.5d0)*amp
       call random_number(tmp) ; vvv(i,j)=2.d0*(tmp-0.5d0)*amp
       enddo
       enddo
                          endif
       enddo
       end subroutine pso_do
!234567890
!      Written by In-Ho Lee, KRISS, August 1, 2018.
       subroutine swarmmove()
       USE modpostest, ONLY : amp,posi,fn
       implicit none
       integer i,j,klo1,klo2
       real*8 test0,tmp,zzr,zzi
       logical loutside

       do klo1=1,nsteps
       do klo2=1,nfactor
       do j=1,npop
       do i=1,ndf
       call random_number(zzr) ; call random_number(zzi)
       vvv(i,j)=omega(j)*vvv(i,j)+cc11(j)*zzr*(pbest(i,j)-posi(i,j))+cc22(j)*zzi*(gbest(i)-posi(i,j)) 
       posi(i,j)=posi(i,j)+vvv(i,j)
       loutside=.false.
       if(abs(posi(i,j)) > abs(amp)) loutside=.true.
       if(loutside)then
       call random_number(tmp) ; posi(i,j)=2.d0*(tmp-0.5d0)*amp
                   endif
       enddo
       test0=fn(posi(1,j))
       if(test0 < epbest(j))then
       epbest(j)=test0 ; pbest(:,j)=posi(:,j)
                            endif
       enddo
       test0=minval(epbest)
       do j=1,npop
       if(abs(test0-epbest(j)) < 1.d-9)then
       gbest(:)=pbest(:,j)
                                       endif
       enddo
       enddo
       enddo
       end subroutine swarmmove
       end module pso
!234567890
       program test_psodriver
       USE pso, ONLY : pso_initial,pso_final,pso_do
       implicit none
       integer ndf0,npop0,mpsoupdate0,nfactor0
       real*8 test

       ndf0=10 ; npop0=20 ; mpsoupdate0=20 ; nfactor0=1
       call pso_initial(ndf0,npop0,mpsoupdate0,nfactor0)
       call pso_do()
       call pso_final(test)
       stop
       end program test_psodriver


  1.058691862857170E-013
   1.00000001481113        1.99999960995437        2.99999999982571     
   4.00000005468264        5.00000014923843        6.00000010936091     
   7.00000121674457        7.99999736929987        9.00000000307783     
   9.99999975757957 


핑백


최근 포토로그