粒子群算法可用于解决强非线性优化问题，原理较为简单(参加：最优化算法之粒子群算法（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, nParticlexPos(:,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 singlecall random_number(rand1); call random_number(rand2(1)); call random_number(rand2)!$OMP end single!$OMP dodo i = 1, nParticle!目标函数objVal(i) = objFun(xPos(:,i))!更新个体最优if(objVal(i)<minErr(i)) thenpBest(:,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 singlegBest = 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
