粒子群算法Fortran代码(OMP并行)

    粒子群算法可用于解决强非线性优化问题，原理较为简单(参加：最优化算法之粒子群算法（PSO）_青萍之末的博客-CSDN博客_粒子群算法)，这里给出Fortran代码实现模块( module POS)。该代码适用于任意参数个数的情况，并采用OpenMP加速计算。

注意：代码求取的是目标函数极小值问题，如需求极大值，需对目标函数objFun取反。

module PSO 
implicit none
contains
!****************************************************************************
! 粒子群优化算法 Particle Swarm Optimization
! 输入:
! objFun     ... 目标函数
! nParticle  ... 粒子数
! maxIts     ... 最大迭代次数
! w          ... 惯性因子(阻尼), 经典值 1.0
! c1, c2     ... 学习因子, 经典值 2.0
! rang       ... 参数范围
! x0         ... 初始值
! 输出:
! x          ... 结果
! err        ... 误差
!****************************************************************************
subroutine particleSwarmOptimization(objFun, nParticle, maxIts, w, c1, c2, rang, x, err, x0)
!$ use omp_lib
implicit none
interface
real function objFun(x) !求极小值
real x(:)
end function
end interface
integer, intent(in):: nParticle, maxIts
real,    intent(in):: w, c1, c2, rang(:,:)
real,   intent(out):: x(:), err
real,    intent(in), optional:: x0(:)
real, allocatable:: xPos(:,:), pBest(:,:), gBest(:), minErr(:)
real, allocatable:: rand1(:), rand2(:), v(:,:), objVal(:)
integer i, its
call RANDOM_SEED()
!参数个数
i = size(x)
allocate(xPos(i,nParticle), v(i,nParticle), pBest(i,nParticle), gBest(i)) !位置、速度, 最佳位置
allocate(rand1(nParticle), rand2(nParticle))   !随机数
allocate(minErr(nParticle), objVal(nParticle))
!初始位置及速度
call random_number(xPos); call random_number(v(1,1)); call random_number(v)
do i = 1, nParticle
  xPos(:,i) = xPos(:,i)*(rang(2,:)-rang(1,:)) + rang(1,:)
  v(:,i) = (v(:,i)-0.5) * (rang(2,:)-rang(1,:))
end do
if(present(x0)) xPos(:,1) = x0(:)

!初始化误差
minErr = huge(1.0)
gBest = xPos(:,1)

!$OMP parallel default(shared) private(fre) num_threads(omp_get_num_procs()-1)
do its = 1, maxIts
  !$OMP single
  call random_number(rand1); call random_number(rand2(1)); call random_number(rand2)
  !$OMP end single
  !$OMP do
  do i = 1, nParticle
    !目标函数
    objVal(i) = objFun(xPos(:,i))
    !更新个体最优
    if(objVal(i)<minErr(i)) then
      pBest(:,i) = xPos(:,i)
      minErr(i) = objVal(i)
    end if
    !更新速度
    v(:,i) = (0.9-0.5*(its-1.0)/max(maxIts-1.0,1.0))*w*v(:,i) + c1*rand1(i)*(pBest(:,i)-xPos(:,i)) &
      + c2*rand2(i)*(gBest(:)-xPos(:,i))
    !限制范围
    v(:,i) = modulo(v(:,i)-rang(1,:),rang(2,:)-rang(1,:)) + rang(1,:)
  end do
  !$OMP end do
  !$OMP single
  gBest = pBest(:,minloc(minErr,1)) !全局最优
  !$OMP end single
end do
!$OMP end parallel
!结果
x = gBest
err = minval(minErr)
deallocate(xPos, pBest, gBest, minErr, rand1, rand2, v, objVal)
end subroutine
end module

测试例子

program Test
use PSO
implicit none
real x(2), rang(2,2), err
!一元函数
rang(:,1) = [0.0,6.0]
call particleSwarmOptimization(fun1, 80, 20, 1.0, 0.5, 1.0, rang(:,1:1), x(1:1), err)
print*,'result: ', x(1:1) !真解 2.0
print*,'error: ', err
!二元函数
rang(:,1) = [-3.0,6.0]
rang(:,2) = [-4.0,3.0]
call particleSwarmOptimization(fun2, 8000, 20, 1.0, 2.0, 2.0, rang, x, err)
print*,'result: ', x!真解 0.5,1.0
print*,'error: ', err
pause
contains
! z= ( x*x + y*y )
function Fun2(x) result(res)
implicit none
real x(:), res
res = (x(1)-0.5)**2 + (x(2)+1.0)**2
end function
! z = (x-2)*(x-2)
function Fun1(x) result(res)
implicit none
real x(:), res
res = (x(1)-2)**2
end function
end program Test