基于改进多目标粒子群算法的环境经济发电调度求解方法

摘要

本发明涉及一种基于改进多目标粒子群算法的环境经济发电调度求解方法，首先建立以燃料费用最低和污染气体排放量最少为调度目标，以发电机出力限制、系统功率平衡及旋转备用为约束的电力系统环境经济调度模型；然后对多目标粒子群算法进行改进，并将改进的多目标粒子群算法应用到电力系统的环境经济调度求解中，同时，根据帕累托占优条件选择个体最优解，采用基于非劣解稀疏度排序的方法选择全局最优解，并根据帕累托最优前沿的斜率特性，提出用斜率法删除远离帕累托最优前沿的解；最后，采用模糊数学法选择帕累托最优前沿的折衷最优解，以辅助决策者选择调度方案。本发明运算简单，运行速度快，且具有较好的稳定性，易于工程实现。

说明书

技术领域

本发明属于电力系统多目标优化调度技术领域，特别涉及一种基于改进多目 标粒子群算法的环境经济发电调度求解方法。

背景技术

电力系统经济调度是在满足电力供需平衡及机组出力上下限的条件下，求解 使系统发电成本或燃料费用最低的调度方案。然而，火电机组在发电过程中不可 避免地向大气排放硫氧化物、氮氧化物及二氧化碳等污染气体。随着人们环保意 识的增强，限制污染气体的排放量也成为电力系统的一个重要调度目标，电力系 统由原来的单目标经济调度转变为多目标环境经济调度。

相比其他减排措施，环境经济调度因其投资少见效快而备受研究者的青睐。 含环境经济调度因素的多目标优化问题，其早期求解方法多采用约束法和权系数 和法。约束法将污染气体排放指标作为经济调度的一个新约束，其求解仍属单目 标优化，但该方法未考虑环境与经济间的折衷关系。权系数和法把多个优化目标 线性组合，从而将多目标优化转化为单目标优化，该方法通过调整各目标的加权 系数获得一组帕累托最优解，但需多次运行算法。另外，该方法缺少选择最优权 系数的基础研究。

环境经济调度是一个高维非线性优化问题，基于梯度等的传统优化方法难以 找到其全局最优解。近年，粒子群算法因此简单、易操作、收敛速度快等特点， 被广泛应用于电力系统的经济调度中，并取得了一定的成效。粒子群算法最早由 美国的Kenndy和Eberhart于1995年提出，用于模拟鸟群的觅食过程，算法通 过粒子间的合作、竞争进行智能指导寻优。目前，已有学者应用粒子群算法来求 解多目标优化问题。将多目标粒子群算法扩展到环境经济调度领域时，算法对多 个目标同时进行优化，但各目标间的矛盾性，使多目标粒子群算法的全局最优解 并不唯一，只能得到一组帕累托最优解。如何合理地选择种群的全局最优解和个 体的最优解成为求解多目标优化问题的关键之一。

发明内容

本发明的目的在于克服现有方法在求解环境经济发电调度方案时系统约束 未能严格满足或得到的帕累托最优前沿（Pareto optimal front，POF）分布不 够广泛均匀的弊端，提出了一种新的基于改进多目标粒子群算法的环境经济发电 调度求解方法，对环境经济发电调度建立了数学模型，给出了使约束严格满足的 处理方法，改进了多目标粒子群个体最优解和全局最优解的选取方法，提出了根 据帕累托最优前沿的斜率特性筛选占优解，实现了帕累托最优前沿分布的广泛性 和均匀性；该方法不仅适合不计网损的环境经济调度，尤其适合采用B系数矩 阵法求解网损的环境经济调度问题。

为实现本发明的目的，本发明提供的技术方案是：

一种基于改进多目标粒子群算法的环境经济发电调度求解方法，包括以下步 骤，

步骤1、建立环境经济调度的数学模型；

环境经济调度的数学模型包括目标函数和约束条件；目标函数包括燃料费用 最低和污染气体排放量最少；约束条件包括机组出力约束、功率平衡约束、旋转 备用约束；

具体形式为：

$\left(\begin{array}{c}F\left(P_G\right)=\underset{i=1}{\overset{N}{\Sigma }}(a_i+b_i{P}_{\mathrm{Gi}}+c_i{P}_{\mathrm{Gi}}^2)\\ E\left(P_G\right)=\underset{i=1}{\overset{N}{\Sigma }}[{10}^{-2}({\alpha }_i+{\beta }_i{P}_{\mathrm{Gi}}+{\gamma }_i{P}_{\mathrm{Gi}}^2)+{\xi }_i\exp \left({\lambda }_i{P}_{\mathrm{Gi}}\right)]\\ s.t.\\ P_{\mathrm{Gi}\min }\le P_{\mathrm{Gi}}\le P_{\mathrm{Gi}\max }\\ \underset{i=1}{\overset{N}{\Sigma }}{P}_{\mathrm{Gi}}-P_D-P_l=0\\ \underset{i=1}{\overset{N}{\Sigma }}{P}_{\mathrm{Gi}\max }-P_D-P_l\ge P_{\mathrm{SR}}\end{array}\right)---\left(1\right)$

其中，P_G为发电机组的有功出力向量；P_Gi为火电机组i的有功出力；N为 火电机组总数；F(P_G)为系统的总燃料成本；E(P_G)为系统的总污染气体排放量； a_i、b_i、c_i为火电机组i的燃料成本系数；α_i、β_i、γ_i、ζ_i、λ_i分别为火电机组 i的污染气体排放量系数；P_Gimin、P_Gimax分别为火电机组i的最小有功出力和最大 有功出力；P_D为系统总负荷需求；P_l为系统总传输网损；P_SR为给定负荷下的系 统备用容量；

步骤2、对多目标粒子群算法进行改进；

对多目标粒子群算法的改进主要体现在约束处理、个体最优解和全局最优解 的选取及帕累托最优前沿的求解上；

约束处理采用惩罚函数和约束修正因子相结合的方法，并通过适应度函数来 体现，适应度函数的具体形式如式（2）所示：

$\min Y=\left(\begin{array}{c}\min [F\left(X\right),E\left(X\right)];\mathrm{ifA}\\ \min [F\left(X\right)+M_{F}P\left(X\right),E\left(X\right)+M_{E}P\left(X\right);\mathrm{if><mo>]</mo>}\end{array}---\left(2\right)\right)$

式中，F(X)、E(X)分别为经济调度和环境调度的原目标函数；M_F、M_E为 惩罚因子；M_FP(X)、M_EP(X)为惩罚项；A表示用约束修正因子法对粒子修正 后，粒子能满足等式约束和不等式约束；B表示修正后粒子仍至少违背一个约束；

多目标粒子群的个体最优解和全局最优解按如下方法选择：

将种群随机产生的初始解作为各粒子的初始个体最优解，之后根据帕累托非 劣解定义来判断本次迭代产生的粒子与上代粒子间的占优关系，若本次迭代产生 的粒子占优上一代的粒子，则将粒子的个体最优解更新为本次迭代产生的粒子； 否则，对粒子的个体最优解不作更新；

粒子群在每次迭代过程中根据占优条件产生全局帕累托最优解集，历史帕累 托最优解集初值等于全局帕累托最优解集初值，之后将每次迭代产生的全局帕累 托最优解集并入历史帕累托最优解集，并根据粒子间的占优关系更新历史帕累托 最优解集；根据稀疏度排序，选择历史帕累托最优解集中稀疏度最大的粒子作为 全局最优解；

改善帕累托最优前沿解的分布：

采用斜率法筛选帕累托最优前沿上的解；对帕累托最优前沿上的解i，设与 其相邻的粒子分别为解i-1和解i+1，按式（3）、式（4）计算粒子之间的斜率：

$k_{i,i+1}=\frac{(f_{2,i}-f_{2,i+1})/(f_{2,\max }-f_{2,\min })}{(f_{1,i}-f_{1,i+1})/(f_{1,\max }-f_{1\min })}---\left(3\right)$

$k_{i-1,i+1}=\frac{(f_{2,i-1}-f_{2,i+1})/(f_{2,\max }-f_{2,\min })}{(f_{1,i-1}-f_{1,i+1})/(f_{1,\max }-f_{1\min })}---\left(4\right)$

式中，k_i,i+1为帕累托最优解i与其相邻解i+1之间的标准化斜率，k_i-1,i+1为与 帕累托最优解i相邻的两个解i-1与i+1之间的标准化斜率；若k_i,i+1＞k_i-1,i+1，则说 明帕累托最优解i离理想帕累托最优前沿较近，保留这样的解；若k_i,i+1≤k_i-1,i+1，

则说明帕累托最优解i偏离理想帕累托最优前沿较远，删除这样的解；

步骤3、根据步骤2中改进的多目标粒子群算法求解环境经济调度，得到环境经 济调度的帕累托最优前沿；

步骤4、在步骤3得到的帕累托最优前沿基础上，采用模糊数学法计算各帕累托 最优解的满意度；

各帕累托最优解在某一维目标函数中的满意度为：

${\mu }_i^m=\left(\begin{array}{cc}1,& f_i^m\le f^{m,\min }\\ \frac{f^{m,\max }-f_i^m}{f^{m,\max }-f^{m,\min }},& f^{m,\min }<f_i^m<f^{m,\max }\\ 0,& f_i^m\ge f^{m,\max }\end{array}\right)---\left(5\right)$

式中，为第i个最优解的第m维目标函数值，m∈{1,2,...,N_obj}，f^m,min、f^m,max分别为第m维目标函数值的最小值和最大值；

各帕累托最优解的满意度为：

${\mu }_i=\frac{\underset{m=1}{\overset{N_{\mathrm{obj}}}{\Sigma }}{\mu }_i^m}{\underset{i=1}{\overset{N_c}{\Sigma }}\underset{m=1}{\overset{N_{\mathrm{obj}}}{\Sigma }}{\mu }_i^m}---\left(6\right)$

式中，N_C为POF上解的个数，选择满意度最大的帕累托最优解作为折衷最 优解，以辅助决策者选择兼顾经济和环境的最优调度方案。

所述步骤3包括以下步骤，

步骤3.1、初始化；

初始化种群的速度和位置，设置最大迭代次数K，给定学习因子c₁、c₂，惯 性权重初值，罚因子值M_E、M_F，根据式（2）中的适应度函数评价初始种群， 初始化各粒子的个体最优解、种群的全局最优解、全局帕累托最优解集及历史帕 累托最优解集；

步骤3.2、根据式（7）计算本次迭代的惯性权重w；

式中，K为最大迭代次数，k为当前迭代次数，w_max＝0.9，w_min＝0.4；

$w=w_{\max }-\frac{w_{\max }-w_{\min }}{K}k---\left(7\right)$

步骤3.3、根据式（8）更新各粒子的速度和位置；

$\left(\begin{array}{c}{v}_{i,j}^{k+1}=w*v_{i,j}^k+c_1{r}_1*(p_{i,j}^k-x_{i,j}^k)+c_2{r}_2*(g_j^k-x_{i,j}^k)\\ x_{i,j}^{k+1}=x_{i,j}^k+v_{i,j}^{k+1}\end{array}\right)j=\mathrm{1,2},...,n---\left(8\right)$

式中，r₁、r₂为[0,1]区间上服从均匀分布的随机数，、分别为粒子i在 第k次迭代中在第j维空间中的飞行速度和位置，为粒子i的个体最优值在第 k次迭代中在第j维空间中的位置分量，为种群的全局最优值在第k次迭代中 在第j维空间中的位置分量；

步骤3.4、修正各粒子的独立机组的有功出力；

步骤3.5、根据式（2）计算各粒子的适应度值；

步骤3.6、寻找本次迭代的全局帕累托最优解集；

步骤3.7、采用斜率法删除历史帕累托最优解集中偏离理想帕累托最优前沿 较远的解；

步骤3.8、在历史帕累托最优解集中，污染气体排放量取最小值和燃料成本 取最小值的两个解除外，对其余解按稀疏度排序，选择稀疏度最大的粒子作为全 局最优解；

步骤3.9、判断历史帕累托最优解集中解的个数是否超过预先设定的容量 N_C，若没有超过，转第十步；若超过，保留污染气体排放量取最小值和燃料成 本取最小值的两个解及稀疏度排序位于前N_C-2位的解，删除其余解，转第十步；

步骤3.10、判断是否达到最大迭代次数K，若无，迭代次数加1，清空全局 帕累托最优解集，转第二步；否则直接输出帕累托最优前沿。

所述步骤3.4的实现过程如下，

独立机组的有功出力用下式计算：

$P_{\mathrm{Gd}}=P_D+P_l-\underset{\underset{j\ne d}{j=1}}{\overset{N}{\Sigma }}{P}_{\mathrm{Gj}}---\left(9\right)$

式（9）中P_Gj为系统中除独立机组外的其他机组出力值，若不考虑式(9)中 的网损P_l，则独立机组的有功出力P_Gd可由已知负荷值和其他机组的出力值直接 求解得到；若考虑网损，则,

$P_l=\left(\begin{array}{cc}{P}_{\mathrm{Gd}}& P_{\mathrm{Ga}}\end{array}\right)\left(\begin{array}{cc}{B}_{\mathrm{dd}}& B_{\mathrm{da}}\\ B_{\mathrm{ad}}& B_{\mathrm{aa}}\end{array}\right)\left(\begin{array}{c}{P}_{\mathrm{Gd}}\\ P_{\mathrm{Ga}}^T\end{array}\right)+\left(\begin{array}{cc}{B}_{0d}& B_{0a}\end{array}\right)\left(\begin{array}{c}{P}_{\mathrm{Gd}}\\ P_{\mathrm{Ga}}^T\end{array}\right)+B_{00}---\left(10\right)$

式中,P_Ga为除去独立机组后的其他(N-1)台发电机组的有功出力向量， P_Ga＝[P_G1,P_G2,...,P_Gd-1,P_Gd+1,...,P_GN]，B_aa为(N-1)×(N-1)矩阵，B_da、B_0a均为(N-1)维 行向量，B_ad为(N-1)维列向量，B_dd、B_0d、B₀₀均为标量；

将式(10)代入式(9)，化简后得到一个关于P_Gd的二次方程：

$a{P}_{\mathrm{Gd}}^2+b{P}_{\mathrm{Gd}}+c=0---\left(11\right)$

式中,

a＝B_dd     (12)

$b=P_{\mathrm{Ga}}{B}_{\mathrm{ad}}+B_{\mathrm{da}}{P}_{\mathrm{Ga}}^T+B_{0d}-1---\left(13\right)$

$c=P_{\mathrm{Ga}}{B}_{\mathrm{aa}}{P}_{\mathrm{Ga}}^T+B_{0a}{P}_{\mathrm{Ga}}^T+B_{00}+P_D-\mathrm{sum}\left(P_{\mathrm{Ga}}\right)---\left(14\right)$

式（14）中，sum(P_Ga)为向量P_Ga中各元素之和，解方程(11)可得到关于P_Gd的 两个解则独立机组的有功出力为：

$P_{\mathrm{Gd}}=\min \{P_{\mathrm{Gd}}^1,P_{\mathrm{Gd}}^2\}---\left(15\right)$

与现有技术相比，本发明具有以下有益效果：

1、本发明运算简单，运行速度快，且具有较好的稳定性，易于工程实现； 其用约束修正因子和罚函数相结合的方法来处理系统约束，并在功率平衡中计及 了网损，采用B系数法求解网损。

2、本发明在多目标粒子群算法中引入了全局帕累托最优解集和历史帕累托 最优解集，重新定义了粒子的个体最优解和种群的全局最优解。

3、本发明能够快速求解出环境经济调度的帕累托最优前沿，且前沿分布广 泛均匀，并能根据模糊数学法为决策者提供折衷最优解。

附图说明：

图1为用斜率法改善帕累托最优前沿解的分布情况的说明图。

图2为本发明的算法流程图。

图3为经济调度收敛曲线。

图4为环境调度收敛曲线。

图5为不同初始解下的经济调度最优值（算例1）。

图6为50次不同初始解下的环境调度最优值(算例1)。

图7为斜率法得到的POF(算例1)。

图8为斜率法得到的POF(算例2)。

图9为稀疏度排序法得到的POF(算例1)。

图10为稀疏度排序法得到的POF(算例2)。

具体的实施方式

下面结合附图对本发明作进一步说明。

本发明包括以下步骤：

步骤1、建立环境经济调度的数学模型

电力系统环境经济调度的数学模型包括目标函数和约束条件两部分。目标函 数包括燃料费用最低和污染气体排放量最少；约束条件包括机组出力约束、功率 平衡约束、旋转备用约束。

目标函数1：燃料费用最低。燃料费用最低属于经济调度，每台发电机组的 燃料成本曲线通常用一个二次函数来表示，系统总燃料成本F(P_G)可表示为：

$F\left(P_G\right)=\underset{i=1}{\overset{N}{\Sigma }}(a_i+b_i{P}_{\mathrm{Gi}}+c_i{P}_{\mathrm{Gi}}^2)---\left(1\right)$

式中：N为火电机组总数；a_i、b_i、c_i为火电机组i的燃料成本系数；P_Gi为 火电机组i的有功出力。P_G为系统火电机组有功出力向量，可表示为：

$P_G=[P_{G1},P_{G2},...,P_{\mathrm{GN}}]---\left(2\right)$

目标函数2：污染气体排放量最少。污染气体排放量最少属于环境调度，火 电机组在发电过程中排放多种污染气体，各污染气体排放量都可与火电机组有功 出力单独建立函数关系，但为方便计算，采用污染气体综合排放模型，则系统总 污染气体排放量可表示为：

$E\left(P_G\right)=\underset{i=1}{\overset{N}{\Sigma }}[{10}^{-2}({\alpha }_i+{\beta }_i{P}_{\mathrm{Gi}}+{\gamma }_i{P}_{\mathrm{Gi}}^2)+{\zeta }_i\exp \left({\lambda }_i{P}_{\mathrm{Gi}}\right)]---\left(3\right)$

式中：α_i、β_i、γ_i、ζ_i、λ_i分别为火电机组i的污染气体排放量系数。

火电机组出力约束:

P_Gimin≤P_Gi≤P_Gimax     (4)

式中:P_Gimin、P_Gimax分别为火电机组i的最小有功出力和最大有功出力。

功率平衡约束:

$\underset{i=1}{\overset{N}{\Sigma }}{P}_{\mathrm{Gi}}-P_D-P_l=0---\left(5\right)$

式中：P_D为系统总负荷需求，P_l为系统总传输网损。网损求解采用B系数法， 其计算公式如下所示：

$P_l=P_{G}B{P}_G^T+B_0{P}_G^T+B_{00}---\left(6\right)$

式中：为P_G的转置，B为N×N维矩阵，B₀为N维行向量，B₀₀为一个标 量。

旋转备用约束:

电力系统存在机组停运及负荷预测误差等问题，为应对这些问题给电力系统 调度带来的影响，在调度中要考虑旋转备用约束。系统旋转备用容量满足如下约 束：

$\underset{i=1}{\overset{N}{\Sigma }}{P}_{\mathrm{Gi}\max }-P_D-P_l\ge P_{\mathrm{SR}}---\left(7\right)$

式中：P_SR为给定负荷下的系统备用容量，一般取系统总负荷的5%。

将上述目标函数及各约束条件整合在一起，即可得到EED （environmental/economic dispatch）问题的数学模型：

$\left(\begin{array}{c}\min Y\left(X\right)=\min [F\left(X\right),E\left(X\right)]\\ s.t.\\ h_i\left(X\right)=0,i=\mathrm{1,2},...,L\\ g_j\left(X\right)\ge 0,j=\mathrm{1,2},...,K\end{array}\right)---\left(8\right)$

式中：X＝P_G为优化问题的解向量；h、g分别为模型的等式约束和不等式 约束；L、K分别为等式约束和不等式约束的总数。

步骤2、对多目标粒子群算法进行改进

选用粒子群算法求解电力系统环境经济调度方案时，首先要将标准的粒子群 算法扩展成多目标粒子群算法，然后对多目标粒子群算法加以改进。对多目标粒 子群算法的改进主要体现在约束处理、个体最优解和全局最优解的选取及帕累托 最优前沿的求解上。

改进约束处理方法：

粒子群算法在求解约束优化问题时，其处理约束的方法主要分为两类：1） 惩罚函数法；2）设计约束修正因子。基于此，本发明将上述两种方法相结合， 提出了一种新的约束处理方法。

首先设计约束修正因子，该约束处理方法仅针对等式约束。

设p_Gi＝[P_G1,P_G2,...,P_GN]代表系统的一个有功调度方案，p_Gi中的各元素分别代 表N台发电机组的有功出力。为满足式(5)中的有功平衡，选N台机组中容量最 大的一台作为独立机组，独立机组的有功出力用下式计算：

$P_{\mathrm{Gd}}=P_D+P_l-\underset{\underset{j\ne d}{j=1}}{\overset{N}{\Sigma }}{P}_{\mathrm{Gj}}---\left(9\right)$

式（9）中P_Gj为系统中除独立机组外的其他机组出力值，若不考虑式(9)中 的网损P_l，则独立机组的有功出力P_Gd可由已知负荷值和其他机组的出力值直接 求解得到；若考虑网损，则对式(6)中P_G、B、B₀进行分割，并将式(6)改写成 如下形式：

$P_l=\left(\begin{array}{cc}{P}_{\mathrm{Gd}}& P_{\mathrm{Ga}}\end{array}\right)\left(\begin{array}{cc}{B}_{\mathrm{dd}}& B_{\mathrm{da}}\\ B_{\mathrm{ad}}& B_{\mathrm{aa}}\end{array}\right)\left(\begin{array}{c}{P}_{\mathrm{Gd}}\\ P_{\mathrm{Ga}}^T\end{array}\right)+\left(\begin{array}{cc}{B}_{0d}& B_{0a}\end{array}\right)\left(\begin{array}{c}{P}_{\mathrm{Gd}}\\ P_{\mathrm{Ga}}^T\end{array}\right)+B_{00}---\left(10\right)$

式中:P_Ga为除去独立机组后的其他(N-1)台发电机组的有功出力向量， P_Ga＝[P_G1,P_G2,...,P_Gd-1,P_Gd+1,...,P_GN]，B_aa为(N-1)×(N-1)矩阵，B_da、B_0a均为(N-1)维 行向量，B_ad为(N-1)维列向量，B_dd、B_0d、B₀₀均为标量。

将式(10)代入式(9)，化简后得到一个关于P_Gd的二次方程：

${\mathrm{aP}}_{\mathrm{Gd}}^2+{\mathrm{bP}}_{\mathrm{Gd}}+c=0---\left(11\right)$

式中：

a＝B_dd

(12)

$b=P_{\mathrm{Ga}}{B}_{\mathrm{ad}}+B_{\mathrm{da}}{P}_{\mathrm{Ga}}^T+B_{0d}-1---\left(13\right)$

$c=P_{\mathrm{Ga}}{B}_{\mathrm{aa}}{P}_{\mathrm{Ga}}^T+B_{0a}{P}_{\mathrm{Ga}}^T+B_{00}+P_D-\mathrm{sum}\left(P_{\mathrm{Ga}}\right)---\left(14\right)$

式（14）中，sum(P_Ga)为向量P_Ga中各元素之和，解方程(11)可得到关于P_Gd的 两个解则独立机组的有功出力为：

$P_{\mathrm{Gd}}=\min \{P_{\mathrm{Gd}}^1,P_{\mathrm{Gd}}^2\}---\left(15\right)$

通过设计约束因子来修正不满足等式约束的粒子，修正后有两种结果，第 一种：等式约束满足，且独立机组出力P_Gd满足其出力约束，即 P_Gd∈[P_Gdmin,P_Gdmax]；第二种：等式约束满足，但独立机组出力P_Gd超出其出力范 围。

当修正后的结果属于第二种情况时，若P_Gd＜P_Gdmin，则令P_Gd＝P_Gdmin；若 P_Gd＞P_Gdmax，则令P_Gd＝P_Gdmax。

当使用约束修正因子法仍不能使粒子满足等式约束或粒子不满足不等式约 束时，则通过罚函数法来处理约束。

对式(8)中的约束优化问题作辅助函数：

Z(X,M_F,M_E)＝[P_F(X,M_F),P_E(X,M_E)]     (16)

其中：

P_F(X,M_F)＝F(X)+M_FP(X)     (17)

P_E(X,M_E)＝E(X)+M_EP(X)     (18)

$P\left(X\right)=\underset{i=1}{\overset{L}{\Sigma }}{\left|h_i\left(X\right)\right|}^2+\underset{j=1}{\overset{K}{\Sigma }}{\left[\min \right\{g_j\left(X\right),0\left\}\right]}^2---\left(19\right)$

P_F、P_E分别是关于燃料成本和污染气体排放量的惩罚函数，M_F、M_E为惩 罚因子，M_FP(X)、M_EP(X)为惩罚项。

取M_F、M_E为适当大小的正数，则式(8)中的约束优化问题转化为求无约束 问题Z(X,M_F,M_E)＝[P_F(X,M_F),P_E(X,M_E)]取最小值的解的问题。M_F、M_E的取值 与目标函数值的数量级相同或相近，其具体取值大小由实验人员根据仿真测试结 果自行调整和选取。

若无约束问题minZ(X,M_F,M_E)的最优解X^*满足等式约束h_i(X)＝0,i＝1,2,...,L 及不等式约束g_j(X)≥0,j＝1,2,...,K，则X^*就是原问题的最优解。

将所提约束处理方法应用到多目标粒子群算法中，即按如下公式计算各粒子 的适应度值：

$\min Y=\left(\begin{array}{cc}\min [F\left(X\right),E\left(X\right)];\mathrm{if}& A\\ \min [F\left(X\right)+M_{F}P\left(X\right),E\left(X\right)+M_{E}P\left(X\right);\mathrm{if}& B]\end{array}\right)---\left(20\right)$

式中：A表示用约束修正因子法对粒子修正后，粒子能满足等式约束和不等 式约束；B表示，修正后粒子仍至少违背一个约束。

改进多目标粒子群个体极值解和全局极值解的选取方法：

多目标优化问题的基本特点之一在于各目标间的矛盾性，即用某种方案改进 某个目标值的同时，该方案可能使另一个目标值变劣。各目标间的矛盾性使多目 标优化问题的解并不唯一，取而代之的所有非劣解的集合，也称为帕累托最优解 集，其定义如下。

定义1：对于求最小值的多目标优化问题，若x₁、x₂都是可行域中的解向量， 当且仅当以下两个条件同时满足时，称x₁占优，或称x₁支配x₂：

$\forall i\in \{\mathrm{1,2},...,N_{\mathrm{obj}}\}:f_i\left(x_1\right)\le f_i\left(x_2\right)---\left(21\right)$

$\exists j\in \{\mathrm{1,2},...,N_{\mathrm{obj}}\}:f_j\left(x_1\right)<f_j\left(x_2\right)---\left(22\right)$

式中，N_obj为目标函数的个数。在目标函数的整个可行域中，若不存在其他 可行解支配x₁，则称x₁为非劣解或帕累托最优解。由所有帕累托最优解组成的集 合称为帕累托最优解集或帕累托最优前沿。

由于多目标粒子群算法并不存在唯一的全局最优解，因此需重新定义多目标 粒子群的个体最优解和全局最优解。本发明提出按如下方法选择多目标粒子群的 个体极值解和全局极值解：

在多目标粒子群算法中设置全局帕累托最优解集和历史帕累托最优解集。全 局帕累托最优解集用于存放当前迭代过程中产生的所有帕累托最优解，可按如下 流程快速求解得到：

假设一个种群中含有m个粒子，每个粒子有N_obj个目标函数值，则通过以下 程序寻找每次迭代产生的全局帕累托最优解集：

1)令i＝1；

2)对所有j＝1,2,...,m且j≠i，用式(21)及式(22)来比较粒子x_i与粒子x_j；

3)若存在j，使粒子x_j支配x_i，则将粒子x_i标记为劣解；

4)若i＞m，转向5）；否则令i＝i+1，转向2）；

5)除去所有被标记的解，剩下的所有解组成本次迭代的全局帕累托最优解 集。

历史帕累托最优解集用来存放整个迭代过程中的帕累托最优解。每次迭代过 程中，将本次迭代产生的全局帕累托最优解集并入历史帕累托最优解集，并根据 式(21)和式(22)中的帕累托占优条件寻找其中的非劣解，删除所有劣解。

个体最优解：将种群随机产生的初始解作为各粒子的初始个体最优解，之后 根据式(21)和式(22)判断本次迭代产生的粒子和上代粒子间的占优关系，若本次 迭代产生的粒子占优上一次迭代的个体最优解，则将粒子的个体最优解更新为本 次迭代产生的粒子；否则，对粒子的个体最优解不作更新。

全局最优解：全局最优解从历史帕累托最优解集中选取。根据解集中各粒子 的稀疏度排序，选择稀疏度最大的粒子作为当前迭代的全局最优解。

改善帕累托最优前沿解的分布：

本发明提出用斜率法来改善帕累托最优前沿解的分布情况。如附图1所示， 以一个双目标优化问题为例，f₁、f₂分别为各目标函数值，假设所得帕累托最 优前沿上有A、B、C、D、E、F6个解。从附图1中可以看出，A、B、D、 E、F5个点均位于理想帕累托最优前沿曲线上（图中虚线所示），C点虽然也 属于帕累托最优解，但其偏离最优前沿，且当C点偏离较远时，有k_CD＜k_BD，k_CD为C、D两点之间的斜率，k_BD为B、D两点之间的斜率。由于f₁、f₂的量纲不 同，因此按下式计算帕累托最优解之间的斜率：

$k_{i,i+1}=\frac{(f_{2,i}-f_{2,i+1})/(f_{2,\max }-f_{2,\min })}{(f_{1,i}-f_{1,i+1})/(f_{1,\mathrm{amx}}-f_{1\min })}---\left(23\right)$

$k_{i-1,i+1}=\frac{(f_{2,i-1}-f_{2,i+1})/(f_{2,\max }-f_{2,\min })}{(f_{1,i-1}-f_{1,i+1})/(f_{1,\max }-f_{1\min })}---\left(24\right)$

式中：k_i,i+1为帕累托最优解i与其相邻解i+1之间的标准化斜率，k_i-1,i+1为与 帕累托最优解i相邻的两个解i-1与i+1之间的标准化斜率。若k_i,i+1＞k_i-1,i+1，则说 明帕累托最优解i离理想帕累托最优前沿较近，保留这样的解；若k_i,i+1≤k_i-1,i+1， 则说明帕累托最优解i偏离理想帕累托最优前沿较远，删除这样的解。

步骤3、用步骤2中改进的多目标粒子群算法求解环境经济调度的帕累托最 优前沿，算法流程如附图2所示，具体步骤如下所示：

步骤3.1：初始化。初始化种群的速度和位置，设置最大迭代次数K，给定 学习因子c₁、c₂，惯性权重初值，罚因子值M_E、M_F，根据式（20）中的适应 度函数评价初始种群，初始化各粒子的个体最优解、种群的全局最优解、全局帕 累托最优解集及历史帕累托最优解集。

步骤3.2：根据式（25）计算本次迭代的惯性权重w，式中，K为最大迭代 次数，k为当前迭代次数，w_max＝0.9，w_min＝0.4。

$w=w_{\max }-\frac{w_{\max }-w_{\min }}{K}k---\left(25\right)$

步骤3.3：根据式（26）更新各粒子的速度和位置。式中，r₁、r₂为[0,1]区 间上服从均匀分布的随机数，分别为粒子i在第k次迭代中在第j维空间 中的飞行速度和位置，为粒子i的个体最优值在第k次迭代中在第j维空间中 的位置分量，为种群的全局最优值在第k次迭代中在第j维空间中的位置分 量。

$\left(\begin{array}{c}{v}_{i,j}^{k+1}=w*v_{i,j}^k+c_1{r}_1*(p_{i,j}^k-x_{i,j}^k)+c_2{r}_2*(g_j^k-x_{i,j}^k)\\ x_{i,j}^{k+1}=x_{i,j}^k+v_{i,j}^{k+1}\end{array}\right)j=\mathrm{1,2},..,n---\left(26\right)$

步骤3.4：根据式（9）至式（15）修正各粒子的独立机组出力。

步骤3.5：根据式（20）计算各粒子的适应度值。

步骤3.6：寻找本次迭代的全局帕累托最优解集。

步骤3.7：用斜率法删除历史帕累托最优解集中偏离理想帕累托最优前沿较 远的解。

步骤3.8：在历史帕累托最优解集中，污染气体排放量取最小值和燃料成本 取最小值的两个解除外，对其余解按稀疏度排序，选择稀疏度最大的粒子作为全 局最优解。

步骤3.9：判断历史帕累托最优解集中解的个数是否超过预先设定的容量 N_C，若没有超过，转第十步；若超过，保留污染气体排放量取最小值和燃料成 本取最小值的两个解及稀疏度排序位于前N_C-2位的解，删除其余解，转第十步。

步骤3.10：判断是否达到最大迭代次数K，若无，迭代次数加1，清空全局 帕累托最优解集，转第二步；否则直接输出帕累托最优前沿。

步骤4、在步骤3得到的帕累托最优前沿基础上，采用模糊数学法计算各帕 累托最优解的满意度。

在得到优化问题的帕累托最优前沿后，采用模糊数学法计算各帕累托最优解 的满意度，以辅助决策者选择折衷最优解。各帕累托最优解在某一维目标函数中 的满意度为：

${\mu }_i^m=\left(\begin{array}{cc}1,& f_i^m\le f^{m,\min }\\ \frac{f^{m,\max }-f_i^m}{f^{m,\max }-f^{m,\min }},& f^{m,\min }<f_i^m<f^{m,\max }\\ 0,& f_i^m\ge f^{m,\max }\end{array}\right)---\left(27\right)$

式中：为第i个最优解的第m维目标函数值，m∈{1,2,...,N_obj}，f^m,min、f^m,max分别为第m维目标函数值的最小值和最大值。

各帕累托最优解的满意度为：

${\mu }_i=\frac{\underset{m=1}{\overset{N_{\mathrm{obj}}}{\Sigma }}{\mu }_i^m}{\underset{i=1}{\overset{N_c}{\Sigma }}\underset{m=1}{\overset{N_{\mathrm{obj}}}{\Sigma }}{\mu }_i^m}---\left(28\right)$

式中：N_C为POF上解的个数，选择满意度最大的帕累托最优解作为折衷最 优解。

实施例

这里应用改进的多目标粒子群算法求解环境经济调度问题，并用IEEE6机 30节点标准系统来验证本发明的有效性。系统总负荷为283.4MW,各发电机组 的出力上下限、燃料成本系数及污染物排放量系数如表1所示。

表1发电机组数据

计算网损用的B、B₀、B₀₀均为系统基准容量100MVA下的标幺值，其值如 下所示：

$B=\left(\begin{array}{cccccc}0.13382& -0.0299& 0.0044& -0.0022& -0.0010& -0.0008\\ -0.0299& 0.0487& -0.0025& 0.0004& 0.0016& 0.0041\\ 0.0044& -0.0025& 0.0182& -0.0070& -0.0066& -0.0066\\ -0.0022& 0.0004& -0.0070& 0.0137& 0.0050& 0.0033\\ -0.0010& 0.0016& -0.0066& 0.0050& 0.0109& 0.0005\\ -0.0008& 0.0041& -0.0066& 0.0033& 0.0005& 0.0244\end{array}\right)$

B₀＝[-0.0107 0.0060 -0.0017 0.0009 0.0002 0.0030]

B₀₀＝9.8573×10^-4

为说明所提算法的有效性，考虑以下两种不同复杂程度的仿真：

算例1：为将所提算法仿真得到的结果与其它文献中的已知结果进行比较， 算例1不考虑功率平衡约束中的网损；

算例2：考虑功率平衡约束中的网损。

为得到帕累托最优前沿的两个边界解，即燃料费用取最小值的解和污染气体 排放量取最小值的解，并验证所提多目标粒子群算法得到的帕累托最优前沿上的 解是否具有良好的分布特性，首先用单目标粒子群算法分别对系统进行经济调度 和环境调度。在单目标粒子群算法中，粒子数设为60；学习因子c₁＝c₂＝2；罚 因子M_F＝10、M_E＝0.0001；最大迭代次数设为100。

算例1和算例2的调度结果如表2、表3所示，从表2和表3可以看出，系 统的等式约束能够严格满足。

表2算例1环境经济调度结果

表3算例2环境经济调度结果

图3和图4分别为经济调度和环境调度的收敛曲线，从图3和图4可以看出， 本发明方法具有较好的收敛性。

为了说明不同初始解对调度结果的影响，让本发明所提算法随机产生50次 初始解，所得经济调度最优值与环境调度最优值的分布分别如图5和图6所示。

从图5和图6可以看出：本发明方法具有良好的全局搜索能力，50次优化均 能找到全局最优值，且不受初始解的影响。值得注意的是，当最大迭代次数减少 时，能够找到全局最优解的次数会相应减少。

在多目标优化中，每代粒子数设为60，最大迭代次数为1000，在算例1和 算例2中，帕累托最优前沿上解的个数都设为30。用所提多目标粒子群算法分 别对算例1和算例2进行优化，所得帕累托最优前沿如图7、图8所示。若仅使 用稀疏度排序法，所得帕累托最优前沿分别如图9和图10所示。

比较图7、图8、图9和图10可以看出，本发明提出的斜率法得到的帕累托 最优前沿要比仅使用稀疏度排序法得到的帕累托最优前沿更光滑均匀。

从图7和图8可以看出，帕累托最优前沿的两个边界点分别对应经济调度和 环境调度的最优值，各边界点的计算结果如表4所示，折衷最优解如表5所示。

表4POF的边界解

表5折衷最优解

机组编号 算例1 算例2 G₁25.1935 25.4726 G₂36.9593 37.6234 G₃53.7415 57.0078 G₄71.0606 68.1824 G₅53.5121 54.6521 G₆42.9330 43.0657 燃料费用/($/h) 608.8184 616.0108 污染气体/(t/h) 0.2015 0.2006

将表4中的结果与表2和表3中的单目标优化结果进行对比，可以看出二者 相差很小，这说明本发明提出的改进的多目标粒子群算法在不同情况下均能找到 较好的边界解，所得帕累托最优前沿分布范围广。