11.5　粒子群算法在科学研究中的应用

例11.1　试用粒子群算法求下列Rastrigin函数的极小值：

解：此函数在n＝2时，有50多个局部极小值，全局极小值为0，位于x_i＝0，其图像见图11.3。

图11.3　Rastrigin函数图像

根据粒子群算法的原理，编制以下函数进行计算：

计算结果如下：自变量值为1.0e－003*（0.0212，-0.1250）时，有极小值2.5530e－006。

例11.2　某水文站有13组水位流量原始观察数据，试用水位流量关系式Q＝aH^b进行拟合。其中Q为流量，H为水位，a、b为常数。

解：在模型参数进行拟合时，常用的方法是最小二乘法。但当观察数据存在极端值时，最小二乘法会产生严重偏差，此时可采用其他的拟合准则，例如残差平方和准则、绝对残差绝对值和准则、相对残差绝对值和准则等。

用例11.1给出的粒子群算法作适当修改，分别采用残差平方和准则、绝对残差绝对值和准则及相对残差绝对值和准则对数据进行拟合，得到结果如下：

绝对残差绝对值准则：a＝4.0287，b＝1.8515。

残差平方和准则：a＝5.3723，b＝1.7325。

相对残差绝对值和准则：a＝4.2657，b＝1.8308。

例11.3　旅行商问题（Traveling Salesman Problem，TSP）一直以来是具有广泛应用价值和重要理论价值的组合优化NP难题之一，常用来验证智能启发式算法的有效性。其形式化定义为：给定n个城市和两两城市之间的距离，求一条访问各城市一次且仅一次的最短路线。现有14城市距离（km）的坐标值，试用粒子群算法求解这个TSP问题。

14个城市距离（km）的坐标为

解：粒子群算法中最关键的是惯性权重的确定，为此提出了各种计算方法。一般来说，当群体的最优适应度值长时间未发生变化（停滞）时，应根据群体早熟收敛程度自适应地调整惯性权重。但若对整个群体采用同样的自适应操作，有可能会破坏优秀粒子群体的性能，从而使算法的性能有所下降。

为了充分发挥自适应调节的效能，应该对处于不同状态的粒子采用不同的惯性权重。根据个体适应度值的不同，可以将群体分为两个子群，分别采用不同的自适应操作，使得群体始终保持惯性权重的多样性。惯性权重较小的粒子用来进行局部寻优，加速算法收敛；惯性权重较大的粒子早期用来进行全局寻优，后期用来跳出局部最优，避免早熟收敛。这样，具有不同惯性权重的粒子各尽其责，全局寻优和局部寻优同时进行，在保证算法能全局收敛和收敛速度之间做了一个很好的折衷。

具体做法是，设粒子群的大小为n，第k次迭代中粒子P_i的适应度值为f_i，最优粒子的适应度值为f_max；粒子群的平均适应度值为。将适应度值高于的适应度值再求平均得到，将适应度值低于的适应度值再求平均得到。

当粒子的适应度值好于时，表示已接近全局最优，此时应被赋予较小的惯性权重，以加速向全局最优收敛，计算公式如下：

其中，w_min为w的最小值，本文取0.10。

当粒子的适应度值差于时，这些粒子为群体中较差的粒子，对其惯性权重的调整采用自适应调整遗传算法控制参数的方法，即如下的计算公式：

其中，k₁为一适当的常数。

根据以上思路进行编程计算，适应度函数为路径的总长，k₁＝0.01。其程序与例11.1中的相似，在此就不再列出。其中一次的计算结果为：1—2—14—3—4—12—6—5—7—13—8—9—11—10—1，其路长为31.8072km，见图11.4。

图11.4　一次的计算结果

例11.4　粒子群算法直观，易于理解，寻优策略简单，调试参数少，收敛速度快且简单易行，因此被广泛应用于求解各种非线性、不可微的复杂优化问题。但是算法本身也有局限性，在优化早期，能迅速向最优值靠拢，但在最优值附近收敛较慢，容易出现所谓的“早熟”，即局部收敛。粒子群算法出现早熟现象的主要原因就是缺乏种群的多样性。

改进的方法是：①引入动态改变惯性权重的概念，使权重随粒子位置的改变而改变，更易摆脱局部极值的干扰；②同时鉴于混沌运动的遍历性特点，将其引入粒子群算法中，简单来说即是混沌初始化和混沌扰动，利用混沌运动的遍历性，产生大量种群，根据粒子间欧式距离，从中提取分布均匀的粒子群初始粒子，使粒子在解空间分布均匀，然后在进化过程中，对粒子最优位置进行局部搜索，提高粒子群算法的开发能力，防止算法早熟，增强算法的全局探测能力。

试用上述思想，对下列函数进行优化。

解：函数的图像见图11.5。

图11.5　函数的图像

（1）动态惯性权重的计算公式如下：

其中，N为粒子数目；为第t次迭代时i粒子的函数值；为t次迭代时最优粒子的函数值。

（2）改进的粒子群算法流程如下：

①混沌初始化。设需要优化的变量为D维。随机产生一个D维向量z₁＝［z₁₁，z₁₂，…，z_1D］，每个分量的范围为［0，1］。然后根据Logistic方程得到M个分量，z₁，z₂，…，z_M：

z_n＋1＝μz_n（1－z_n），　　n＝0，1，2，…；　0＜z_n＜1；　μ∈［0，4］

并根据下列公式将混沌区间映射到变量的取值范围：x_ij＝a_j＋（b_j－a_j）z_ij，其中b_j、a_j分别为优化变量的上下限，i＝1，2，…，M；j＝1，2，…，D）。根据目标函数计算每个粒子的适应度值，从M个初始粒子群中选择性能较好的N个作为初始解，随机产生粒子的速度。

②设置粒子的初始个体极值和全局极值。定义各粒子的当前位置为个体极值P_i，根据目标函数计算各个体极值P_i对应的适应度值，选其最好值的粒子位置定义为全局极值P_g。

③根据速度和位置的更新公式更新粒子的飞行速度和位置。

④对最优位置P_g进行混沌优化。先将最优位置映射到Logistic方程的定义域［0，1］，即使用公式z_g＝（P_g－a_i）／（b_i－a_i），再根据Logistic方程进行迭代产生m个混沌变量序列，最后把产生的混沌变量序列映射返回到优化变量的取值区间，得到m个粒子，计算每个粒子的适应度值，得到最优解p′。

⑤用p′代替当前群体中任意粒子的位置。

⑥返回步骤③，直到满足粒子群终止条件则停止计算，输出计算结果。

根据改进的粒子算法流程，编制相应程序进行计算。限于篇幅，在此只列出混沌化的程序，其余与基本的粒子群算法相差不大。

经试验，迭代计算20次左右就可以得到较为满意的结果：

a＝1.0e－003*（0.0767－0.4004），b＝2.6414e－005

例11.5　求下列函数的最小值：

解：此函数的图像见图11.6。

图11.6　函数的图像

本例是带约束条件的优化问题。对于求解这类问题的核心就是如何处理约束条件。传统的方法有梯度映射法、梯度下降法、惩罚函数法和障碍函数法等。

常用的是惩罚函数法，即先是不考虑约束条件地去产生潜在解，然后通过降低评价函数的“好坏度”对其进行惩罚。惩罚方法有两种：一种是通过引入对违反约束条件的惩罚将约束问题转换为非约束问题，这些惩罚包含在评价函数里，且惩罚函数可以是多样的；另一种惩罚方法是从群体中消去不可行解，即最严重的死亡惩罚。

本文采用如下的惩罚函数作为PSO算法的适应度函数：

其中，∣g_j（x）∣表示当g_j（x）≥0时（不满足约束条件）取∣g_j（x）∣的值，当g_j（x）＜0时（满足约束条件）取0；f_max是群体中最差可行粒子的适应值，群体中没有可行粒子时为0。进化群体中不可行粒子在的选择压力下逐渐向可行域靠近，进入可行域后在f（x）作用下接近最优解。计算结果如下：x₁＝0.5000，x₂＝0.7071时，函数值为0.2500。限于篇幅，下面只列出计算适应度值的程序。

例11.6　某镇是一个重要的工业化城镇园区，园区内工厂企业、学校、科研单位较多，且大多涉及化学品或化学有害物；并且紧挨高速公路，运输化学品（化学危险品）、油料和天然气等的车辆较多。该镇人口较多且较为集中，是某城市重要的二级水源地；毗邻重要的历史文物旅游胜地，周围有大片的农作物作业区。一旦有化学物环境污染突发事件，带来的严重后果可想而知。因此，计划在此建立突发应急服务点，应对可能发生的突发事故，保障人民生命财产安全，保护生态环境良好。图11.7所示为该镇的布局坐标图，图中“·”代表化学品环境污染应急救援点。

图11.7　该镇的布局坐标图

根据地图比例尺，将实际距离缩小，得到各应急点的坐标（坐标数字单位全为mm）如表11.1所列，其比例尺为1（mm）∶10（m）。

表11.1　各应急点的坐标值

解：首先对各应急点进行分析，采用AHP法确定其在救援过程的相对重要性，即确定各应急点的权重。通过计算可以得到表11.2所列的权重。

表11.2　各应急点的应急权重

考虑到平面选址，用二维坐标表示应急救援点及应急服务点。应急服务点坐标用（x，y）表示，其中x为直角坐标系中的横坐标，y为纵坐标；各应急点v_i用（x_i，y_i）表示，则模型为

其中，λ是常数，即应急服务点至图中最远的点距离。考虑实际意义，以110出警为依据（目前国内尚无关于化学环境污染应急时间的规定）：“110报警服务接到群众报警，处警人员在市区，必须5分钟内到达现场；在郊区，必须10分钟内到达现场。”结合本文实例，对于小规模的城镇，以处突时车速60km/h计（扣除交通、车辆本身等因素），设应急范围λ＝2.5km。按照比例尺，在图中λ＝250mm。

本例对约束条件采取如下惩罚函数

结合本例实际要求，最终可得罚函数为

根据以上模型，利用粒子群算法进行优化计算，得到如下的结果：应急服务点的坐标为（115.5419，52.293），此时距离为78.8511。

例11.7　对某地区井田的煤层地质条件进行分析，得到表11.3所列的数据，试利用粒子群算法对其进行分类。

表11.3　各煤层块段的特性指标值

解：在基于粒子群算法的聚类分析中，每个粒子代表K个类的中心点，每个粒子的结构可以表示为z_i＝（c_i1，c_i2，…，c_iK），其中c_ij代表第i个粒子的第j类中心点坐标向量，它具有n维（即样品特征数）。粒子群由许多候选分类方案构成。对分类方案优劣进行评价是应用优化算法进行聚类的关键。一般可以采用如下的适应度函数：

这一适应度函数代表了所有类的类内距离和。虽然聚类的结果与类内距离和类间距离具有一定的关系，但这种关系是不确定的，导致应用这个适应度函数来对聚类结果进行评价策略不十分充分。

本例对样品的聚类由下列适应度函数f决定：

其中，w₁和w₂为适当的正常数；为聚合类c_ij中元素的个数，代表z_i对应分类的最大的类内平均距离，代表z_i对应分类的最小的类间距离。这样，通过搜索f的最小值，可以使分类方案同时满足类内距离小和类间距离大。调整w₁和w₂可以方便地给出不同的优先搜索策略。

利用基本的粒子群算法对表中的数据行聚类分析，其中类别数为3，粒子个数取为30。w₁和w₂分别取0.2和0.6。计算结果如下：第2、3、18个样品可以分成一类，其余各样品为另一类。